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ABSTRACT 


A numerical method is presented to implement structural design 
sensitivity analysis using the versatility and convenience of existing 
finite element structural analysis program and the theoretical 
foundation in structural design sensitivity analysis. Conventional 
design variables, such as thickness and cross-sectional areas, are 
considered. Structural performance functionals considered include 
compliance, displacement, and stress. It is shown that calculations can 
be carried out outside existing finite element codes, using 
postprocessing data only. That is, design sensitivity analysis software 
does not have to be imbedded in an existing finite element code. 

The finite element structural analysis program used in the 
implementation presented is IFAD. Feasibility of the method is shown 
through analysis of several problems, including built-up structures. 
Accurate design sensitivity results are obtained without the uncertainty 
of numerical accuracy associated with selection of a finite difference 


perturbation. 



TABLE OF CONTENTS 


Page 

LIST OF TABLES % iv 

LIST OF FIGURES * . v 

LIST OF SYMBOLS vi 

CHAPTER 

I. INTRODUCTION 1 

1 .1 Purpose 1 

1.2 Adjoint Variable Method 4 

1.3 Adjoint Variable Method Results 4 

II. DESIGN SENSITIVITY ANALYSIS METHOD 6 

2.1 Calculation Procedure for Structural 

Components 6 

2.1.1 Membranes 8 

2.1.2 Bending of Beams 15 

2.1.3 Bending of Plates 20 

2.2 Calculation Procedure for a Built-Up Structure . * 26 

III. NUMERICAL EXAMPLES ......... 36 

3.1 Membranes 36 

3 .2 Bending of Beams 41 

3.3 Bending of Plates 45 

3.4 Built-up Structure 50 

IV. CONCLUSIONS 55 

REFERENCES . 56 

APPENDIX 

DESIGN SENSITIVITY IFAD PROGRAM 57 


iii 



LIST OF TABLES 


Table Page 

1 . Membrane Design Sensitivity Check for Compliance 38 

2. Design Sensitivity Check for Displacement .... 39 

3. Design Sensitivity Check for Stress 40 

4. Beam Design Sensitivity Check for Compliance ... 42 

5. Beam Design Sensitivity Check for Displacement 44 

6. Beam Design Sensitivity Check for Stress 44 

7. Plate Design Sensitivity Check for Compliance 47 

8. Plate Design Sensitivity Check for Displacement ..•••• 47 

9. Plate Design Sensitivity Check for von Mises' Stress .... 49 

10. Built-up Structure Design Sensitivity Check for 
Compliance with 6h = O.Olh, 6b = 0.01b, and 

6t = O.Olt 52 

11. Built-up Structure Design Sensitivity Check for 
Displacement with 6h = O.Olh, 6b = 0.01b, 

and 6t = O.Olt 52 

12. Built-up Structure Design Sensitivity Check for 
von Mises' Stress with 6t = O.OOlt, 

6b * 0.001b, and 6h = O.OOlh 54 


LIST OF FIGURES 


Figure Page 

1. Flow Chart of Feasibility Check Procedure 3 

2. Flow Chart of Design Sensitivity Calculation Procedure • % . 7 

3. Clamped Elastic Solid of Variable Thickness h(x) 8 

4. Cantilever Beam with Variable Width and Height 15 

5. Clamped Plate of Variable Thickness t(x) 20 

6. Roof Structure Example 28 

7. Plane Elastic Solid Finite Element Model 37 

8. Cantilever Beam Finite Element Model 42 

9. Bending Plate Finite Element Model 46 

10. Built-up Structure Finite Element Model 51 


v 



LIST OF SYMBOLS 


u 

Z 




6u 

h 

t 

b 

E 

v 

D(u) 

J 

w,w £ 

G 

6 


e,e 


ij 


a, a 


A, A 


ij 

(i) 


Q 


Bilinear form which is dependent on u 
Load linear form which is dependent on u 
Design vector 

Design sensitivity vector; Gauss quadrature counter; 
Local coordinate system of a beam 

Various constraint fuctionals 

Design sensitivites of various constraint functionals 

Perturbed design vector 

Thickness of membrane; heigth of beam 

Thickness of bending plate 

Width of beam 

Young's modulus 

Poisson's ratio 

Flexural rigidity of plate 

Torsional rigidity of beam; determinant of the Jacobian 

Gauss quadrature weight factor 

Shear modulus 

Kroneker Delta 

Strain 

Stress 

Adjoint variable vector 
Domain of the system considered 


vi 



r 

l 


T.T 1 


6h 

<5h 

St 

p 

3 

[ 1 
d 


Z • j | z % z 
ii* xx* xx. 

i 

z xx* z ii 

k 

Y 

X ,X 

XX ’ XX, 

l 

L 

A 

<fr 


Boundary of the system considered 

Summation 

Displacement 

Applied body force 

Applied traction 

Virtual displacement 

Perturbation of h 

Perturbation of b 

Perturbation of t 

Characteristic function 

Partial derivative 

Matrix 

Virtual displacement of an element 

Beam curvature due to z 

Beam curvature due to z 

Element counter 

Material density 

Beam curavature due to X 

Length of a beam element 

Area of a triangle 

General function 


vii 



1 


CHAPTER I 
INTRODUCTION 


1.1 Purpose 

To date there exists a wide variety of finite element structural 
analysis programs that are used as reliable tools for structural 
analysis. They give the designer pertinent information such as 
stresses, strains and displacements of the mechanical system being 
modeled. However, if this information reveals that the mechanical system 
does not meet specified constraint requirements, the designer must make 
intuitive guesses as to how to improve the design. If the mechanical 
system is complex, it becomes very difficult to decide what step must be 
taken to improve the design. There is however, substantial literature 
on the theory of design sensitivity analysis, which predicts the effect 
that structural design changes have on the performance of a mechanical 
system. Use of this technique has been primarily confined to papers in 
structural optimization literature. 

The purpose of this work is to develop and implement structural 
design sensitivity analysis using the adjoint variable method that takes 
advantage of the versatility and convenience of an existing finite 
element structural analysis program and the theoretical foundation in 
structural design sensitivity analysis that is found in Ref. 1. The 
finite element program that will be used is IFAD [3]. It is developed 
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by Applicon Inc. and has been provided to the Center for Computer Aided 
Design for the use in this study. 

In order to check the feasibility of using the design sensitivity 
analysis technique with IFAD, an approximation of the differential tp* of 
a structural performance measure ij> is made using the finite difference 
method. An appropriate design perturbation 6u must be selected in order 
to insure accuracy of the perturbation A\p of the constraint 
functional. If 6u is too small, AiJ> = u + 6u) - i|»(u) may be inaccurate 
due to loss of significant digits in the difference. On the other hand, 
if 6u is too large, Ai|i will be influenced by nonlinearities and the 
differential approximation will be inaccurate. The feasibility check 
procedure is outlined in the flow chart of Fig. 1. Details of the 
calculations of the constraint functionals, the adjoint loads, and the 
design sensitivity vectors for each constraint functional are described 
in Chapter II, for different types of finite elements. The design 
sensitivity i|>' of the constraint functional is the scalar product of the 
design sensitivity vector £ and the design variable perturbation 
vector 5u. If the design variable is constant throughout the finite 
element model of the mechanical system, this becomes a scalar 
multiplication. If the design sensitivity is an accurate prediction of 
the performance of the mechanical system due to a design change, it 
should be equivalent to the difference of the constraint functionals of 
the two finite element models, the original and the perturbed model. 
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Figure 1. Flow Chart of Feasibility Check Procedure 
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1.2 Adjoint Variable Method 

A number of methods could be used to implement structural design 
sensitivity analysis with an existing finite element code, but the most 
powerful is the adjoint variable method. This method can be implemented 
outside of an existing finite element code, using only postprocessing 
data. This is convienent, because the source code for most finite 
element programs is not readily available. If the code is available, 
less programming is involved. The same subroutines for the element 
shape functions used in the finite element model can be used in the 
design sensitivity analysis, since this method is dependent on element 
type. Generality is another factor that adds flexibility to the adjoint 
variable method. The code can be written to include basic design 
variables, constraints and loading conditions. This enables the 
designer to choose what design variables to modify to give the best 
design improvement. 

The adjoint variable method can easily be used for complex 
mechanical systems that have more than one structural component. The 
details of this procedure are discussed in Section 2.2. Design 
sensitivity of a built-up structure is formed by combining the design 
sensitivities of each structural component. The only precaution that is 
necessary is in making sure that the interaction between the components 
is taken into account. 


1.3 Adjoint Variable Method Results 
The design sensitivity vector is the derivative of the constraint 
functional with respect to the design variables. It has the same number 
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of components as there are elements in the finite element model. The 
magnitude of each component reflects how sensitive the element is to a 
change in design relative to the constraint functional. If the vector 
component is negative, the corresponding design variable should be 
decreased to increase <|>. Likewise, if the vector component is positive, 
the design variable should be increased to increase i|». In addition, if 
the magnitude of the vector component is large, then the corresponding 
design variable will have a more substantial effect on design 
improvement . 

When a designer uses a finite element structural analysis in design 
of a mechanical system, it is most likely that a number of program runs 
are necessary before a substantially improved design is obtained. With 
the aid of a design sensitivity vector, the designer will know what 
direction to take to improve the design most efficiently. 
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CHAPTER II 

DESIGN SENSITIVITY ANALYSIS METHOD 

2.1 Calculation Procedure for Structural Components 
To implement the adjoint variable technique of design sensitivity 
analysis, the adjoint load for each constraint functional must be 
calculated. This procedure is developed in Ref. 1 using compliance, 
displacement, stress and natural frequency as constraint functionals. 

For the compliance functional, the adjoint equation is the same as the 
state equation. In this special case the adjoint system does not need 
to be solved. For the displacement functional the adjoint load is a 
unit point load acting at the point where the displacement constraint is 
imposed. To calculate the adjoint response it is necessary only to 
restart the finite element analysis with unit loads applied at varying 
points along the structure. For the stress functional the shape 
function of the structural component used in the finite element analysis 
must be known. This shape function is used to calculate the adjoint 
load for a stress constraint of a specific element of the structure. 

From this point the procedure is similar to the displacement functional, 
in that a restart of the finite element analysis must be completed using 
the adjoint loads of elements as other load cases. 

The flow chart of Fig. 2 shows the overall process. This procedure 
is implemented after the structural response of the finite element 




Figure 2. Flow Chart of Design Sensitivity Calculation 
Procedure 
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model due to the original load has been solved. The original structural 
response plus the adjoint response for each constraint is then utilized 
to calculate the design sensitivity vectors. 

The following sections give detailed explanations of the 
calculation procedures and equations necessary for analyzing membranes, 
bending beams and bending plates. 


2.1.1 Membranes 

Consider a variable thickness thin elastic clamped solid, as shown 
in Fig. 3. The design variable is taken as the variable thickness 
u = h(x) of the plate. 




Figure 3. Clamped Elastic Solid of Variable Thickness h(x) 


The energy bilinear form and the load linear form of the plane 
elasticity problem are given as [1] 


a (z,z) = // h(x) l o i ^(z)e 1 ^(z)dJ2 
U i,j=l 


( 1 ) 


and 
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2 2 

A (z,i) = fJ Q h(x)[ l F i z i ]d« + / [ l T i z i ]dr (2) 

u “ i=l i=l 

where z = [z*,z^]^ is the displacement, F = [F*,F^]^ is the applied body 

i9t ii xi — 

force, T = [T i ,T z ] i is the traction, and a (z) and e (z) are the 

stress and strain fields associated with the displacement z and the 

virtual displacement z respectively. The state equation is given as [1] 


a u ( 2 * z ) = 


(3) 


for all kinematically admissible virtual displacement z. 

First consider the functional representing the compliance of the 
structure as 


*, = /Jo Mx) [ l F^dft + Jr [ l T i z 1 ]dr (4) 

1 i=l i=l 

The first variation of Eq. (4) is 

♦i - j/o [ l F^Sh dfi+ // h[ \ pV'jdB 
1 “ i=l i=l 

(5) 


+ f T [ I TV’jdr 

i=l 

In order to eliminate the dependence on the state variable in 
Eq. (5), it is necessary to define the adjoint equation as [1] 


a u (X,A) 


//q *[ I 

i=l 


F 1 ^ 1 ]d « + 


/ [ l T^jdr 

i=l 


( 6 ) 


for all kinematically admissible virtual displacement X. Since Eq. (6) 
is identical to Eq. (4) if X = z and X = z, the adjoint equation does 



not need to be solved. Using the adjoint variable method of design 
sensitivity analysis gives [1] 
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r, -Hail rVjSh dn+ I Jo [ l - l o lj (z)e i:i (X)]5h d fl 
1 “ i=l i=l i,j=l 

2 2 

= //o [2 l F i z i - I a ij (z) e ij (z)]6h dQ (7) 

i=l i,j=l 

since z = X for the compliance functional* 

To numerically integrate Eqs. (4) and (7), a two-point Gauss 
quadrature formula is used. The equations become 


N 2 2 2 

I l\[ l l F ( z tl M ( I + l T 1 * 1 ) 

k=l fc=l i=l i=l 


( 8 ) 


and 



N 2 2 


i { n i 2f 

k=l Jt=l i=l 


i i 
l Z % 


2 

- I a^ 3 (z)e^ j (z)W )l ]}j6h k 
i,j=l 


(9) 


respectively, where J is the Jacobian, N is the total number of 
elements, subscript i is the counter for the number of Gauss points, 
subscript k is the counter for the element number, W is the weighting 
constant for the £th Gauss point, and supercript i is the direction of 
the force and the displacement* 

Next consider the functional representing the displacement z at a 
discrete point x as 

\|>2 = z(x) = //g 6 (x - x)z(x)d!J 


( 10 ) 
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where 6(x) is the Dirac measure in the plane, acting at the origin. The 
first variation of Eq. (10) is 

i |»2 = J7q $(x - x)z ' (x)dfl (11) 

The adjoint equation in this case is [1] 

a u (*.X) s //^ 5(x - x)X(x)dfl (12) 

for all kinematically admissible virtual displacement X. This equation 

( 2 ) ( 2 ) 

has a unique solution X , where X is the plate displacement due to 

a 

a unit point load acting at a point x. Using the adjoint variable 
method of design sensitivity analysis gives 

*2 “ //« t I fV 2) - l a ij (z)E ij (X (2) ) ]6h d a (13) 

i=l i,j = l 

( 2 ) 

where X is the solution of Eq. (12). 

For this constraint two equations must be solved. The adjoint load 
of Eq. (12) is just a unit load applied at a discrete point in the 
finite element model. All that is necessary is a restart of the model 
so that load cases of applied unit loads at various nodal points can be 
analyzed. The resulting strains due to the adjoint load are then used 
in calculating Eq. (13). Note that for each displacement constraint 
there is a different adjoint load. 

With numerical techniques applied as in the compliance constraint 
case, Eq. (13) becomes 
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N 2 2 


mu* 

k=l A =1 i=l 


w 1 

£ £ 


- \ o lj ( Z )e ij (X (2) )]w }J 6h 

i.J-l * * ^ k 


(14) 

Finally consider the general functional representing a locally 
averaged stress on an element as 


+* 3 JJ fl g(a(z))m n dft 


(15) 


where m„ is a characteristic function defined on a finite element ft as 
P p 


m = 
P 


T^r • 

P 


x e ft 


x ^ ft 


(16) 


The first variation of Eq. (15) is 


♦3 ‘ l! a I ? , TjT ,,1J<2 ' ) K d!! 


(17) 


i,j=l 3 a 


Replacing the variation in state z’ by a virtual displacement X, the 
adjoint equation is obtained as [1] 


a <*.*) = IJq [ I “%■ o ij (X)]m dft ( 18 ) 

u “ i,j = l 3a 1J p 

for all kinematically admissible virtual displacements X. Eq. (18) has 

(3) 

a unique solution for a displacement field X . Using the adjoint 
variable method of design sensitivity analysis gives 

*3 - JJ Q [ l fV 3) - \ ff ij (z)e lj U (3) )]6h dft (19) 
i“l i,j=l 

(3) 

where X is the solution of Eq. (18). 
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With numerical techniques applied as in the displacement constraint 
case, Eqs. (15), (18), and (19) become 


2 

^3=1 [g (**(*) )« W J J 
J £=1 * P * 


( 20 ) 


// fl t I o 11 (»)].- da 

i,j=l 9 o 1J p 


-Ml « t j 3 

A=1 i,j = l 9 cm p * 


( 21 ) 


and 


N 2 2 


. i 2 


*3 


I U 1 ^ 

( 22 ) 


k=l i=l i=l 


respectively. The numerical calculation of the adjoint load is 
considerably more difficult in the stress functional case. The shape 
function of the element must be known so that a*^(X) can be calculated 
using the finite element technique [2] 

a = [E][B]d = [E] £ (23) 

where [E ] is the elasticity matrix, [B] is the strain-displacement 
matrix, which relates to the element shape function, d is the 
displacement vector, and e is the strain vector. 

For a plane stress problem. 



1 

v 


v 0 

1 0 

1-v 

2 


( 24 ) 


0 


0 
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The adjoint load becomes 

J/ 8 [ l 4 « (25) 

U i,j=l 3o J P 

where F is the adjoint equivalent nodal force. Since ra p is a 
characteristic function on the finite element Q p , F acts only on the 
nodal points of element 8 p . 

After the adjoint load is calculated for various elements, a 
restart of the finite element model for each load case corresponding to 
each element adjoint load is made. The strains resulting from these 
adjoint loads are then used in calculating Eq. (22) for the sensitivity 
of each functional. 

When principal stress is selected as the functional, 


= ( ° U + ° 22)/2 + T max (26) 

t = {[(o 11 - o 22 /2] 2 + (O 12 ) 2 } 172 (27) 

max 1 L J J 

and 

v 3 ” 11 -W < o11 - ° 22)/ w 

Sgj/So 22 “ 7 - J (o 11 - < ’ 22 )/ T ma x < 28 > 

. 12 12 , 

3g./3a = a /x 

& 1 max 

and when von Mises ' stress is selected as the functional, 

r , 11.2 1 1 22 . m 22.2 . 12.2 -,1/2 

gj = [(o ) - o O + (o ) + 3 (a ) J 


( 29 ) 
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and 

8g 2 /3<J ll = I (2 ° U “ ° 22)/g 2 

3g 2 /3a 22 = y (2a 22 - a U )/g 2 (30) 

3g 2 /3 a 12 = 3a 12 /g 2 

2.1.2 Bending of Beams 

Consider a cantilever beam with variable width and height and self 
weight, as shown in Fig. 4. The width and height are the design 
variables, u = [b(x),h(x)] . 

The energy bilinear form and the load linear form of the beam are 



Figure 4. Cantilever Beam with Variable Width and Height 
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where y is the weight density of the beam material, F is the distributed 
load, E is the modulus of elasticity of the beam material, bh^/12 is the 
moment of inertia, z is the virtual displacement, z xx is the beam curva- 
ture, and z xx is the beam curvature due to the virtual displacement z. 

The negative sign in the load linear equation is due to the fact 
that the load is applied in the -z direction. 

The state equation is [1] 

a u (z,z) - £ u (z) (33) 

for all kinematically admissible virtual displacements z. 

First consider the functional representing the compliance of the 
structure as 

* - /q (F + ybh)z dx (34) 

The first variation of Eq. (34) is 

* - /q (F + ybh)z' d x - /q - hyz dx 6b - /q byz dx 6h 

(35) 

To replace the variation in state z' by a virtual displacement X, the 
adjoint equation is defined as [1] 

a u (X,X) = - /q (F + ybh)X dx (36) 

for all kinematically admissible virtual displacements X. Since 
Eq. (36) is identical to Eq. (34), if X = z the adjoint equation does 
not need to be solved. Using the adjoint variable method of design 
sensitivity analysis gives 
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- /g [-2yhz - (Eh 3 /l2)( Zxx ) 2 ]dx 


5b 


(37) 


+ Jjj [-2ybz - (3Ebli /12)(z )' ]dx 5h 


To numerically integrate Eqs. (34) and (37), a three-point Gauss 
quadrature formula is used. These equations become 


♦* * A KV T t+ y \\ ,z i v i 


k-1 4-1 


and 


(38) 


N 3 o o 

K = I { I [- 2 Vh k z,- (Eh /12)(z ) ]W }J 5b 

4 k=l 1=1 K x R xx £ A K 

3 3 

+ l { I [-2Yb k z A - (3E b k h 2 /12)( Zxx ) 2 W }J 5h k 
k=l i=*l j l 


(39) 


where N is the total number of elements, SL is the Gauss point counter, 

W is the weighting constant for the Ath Gauss point and J is the 
Jacobian. The beam curvature z xx is calculated using a cubic polynomial 
for the standard beam shape functions. Because the load is in the -z 
direction, it is necessary to change the sign of the local element y- 
rotation 9 . 

y 

Next consider the functional representing the displacement z at a 
discrete point x as 


= z(x) = /q 5(x - x)z(x)dx 


(40) 


where 6(x) is the Dirac measure at zero. The first variation of 


Eq. (40) is 
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^5 = f L 0 ~ x)z ' (x)dx 


The adjoint equation is defined as [1] 


a u (A,X) = J 0 <5(x - x)X(x)dx 


for all kinematically admissible virtual displacements X. Equation (42) 
has a unique solution X^\ where X^^ is the beam displacement due to a 

A 

unit point load acting at a point x. Using the adjoint variable method 
of design sensitivity analysis gives 


*5 - Jo [-hYX (5) - (Eh 3 /12)z xx X ( ^]dx 6b 


+ fl [-byX° ; - (3Ebh /I2)z X° ; ]dx <5h 

J 0 1 xx xx J 


As in the membrane displacement constraint case, only Eq.(43) needs 
to be solved numerically. Using the three-point Gauss quadrature 
technique, Eq. (43) becomes 


+5- I { I [-\YX (5) - (Eh 3 /12)z X< 3) ]w }J 6b 

3 k=l i=l K * K xx Jt * X 


- X l**^^™*# 6h k 


Finally consider the functional representing the allowable 
stresses in the beam as 


“ 7T hEz ra dx 
6 J 0 2 xx p 


( 45 ) 
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where h/2 is the half-depth of the beam, and ra p is a characteristic 
function defined on a finite element dx p as 



x e dx 

P 


x i dx 

P 


The first variation of Eq. (45) is 


(46) 


*6 ‘ Jo I* 2 hE Z kV 2 E z x*"p 5h ] d,t W) 

Replacing the variation in state z' by a virtual displacement X, the 
adjoint equation is defined as 

a (X,X) = ~ /n T hEX m dx (48) 

u '0 2 xx p 

for all kinematically admissible virtual displacements X. Equation (48) 
has a unique solution for a displacement field X ^ . Using the adjoint 
variable method of design sensitivity analysis gives 

% - /o [-hv* (6) - (Eh 3 /12)z xx i ( ^]dx «b 

(49) 

+ /i [“ T Ez m ~ byX (6) - (3Ebh 2 /12)z X (6) ldx Sh 

J 0 1 2 xx p xx xx J 

where X^^ is the solution of Eq. (48). 

With the three-point Gauss quadrature numerical integration 
technique, the integrals in Eqs. (45), (48), and (49) become 


* 


6 


N 3 

HI (- 


k=l A= 1 


1 u „ 

■=- n. E z m 

2 k xx z p 



( 50 ) 
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and 


fb v E X m dx 
J 0 2 xx p 


? (I l- T e[b] » Jw )j d 

k=i z * p * 


( 51 ) 


K " i Jj f " ( Eh k/ 12 ) z xx /i x ; J w i } J - k 




,(6) 

,1 xx . 
i SL 


5b, 


v! u [-7®^" n -v^ 6) 


k=l Jt=l 


it ui u, l /v - 

xx^ p k i, 


(52) 


- (3e 51^ 

where [B] is the strain-displacement matrix, which is the second 
derivative of the shape functions. 


2.1.3 Bending of Plates 

Consider the clamped plate in Fig. 5 of variable thickness 
u = t(x), with a distributed load f(x) that consists of an externally 
applied pressure F(x) and self weight, given by [1] 

f(x) = F(x) + yt(x) (53) 

where Y is the weight density of the material. 



Figure 5 


Clamped plate of variable thickness t(x) 
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For this design dependent loading, the energy bilinear form and the 
load linear form for the plate are given as [ 1 ] 

a 

a u ( z,z) = // fi D(u) [z 11 z 11 + z 2 2 z 22 + V ( Z 22 Z 11 + Z 11 Z 22^ 

+ 2(1 - v)z 12 z 12 ]dSl (54) 

and 

VS) = JJq [F + Yt]z d a (55) 

where D(u) = Et7[12(l-v‘')] is the flexural rigidity, E is Young's 
modulus, v is Poisson's ratio, F is the externally applied pressure, 
and Y is the material density. The state equation is [1] 

a u (z,z) = V z ) (56) 

for all kinematically admissible virtual displacements z. 

First consider the functional representing the compliance of the 
structure as 

* 7 - JJ a <F + Yt)z dn (57) 

The first variation of Eq. (57) is 

♦7 - JJ Q [ (F + Yt)z' + Yz St]dG (58) 

The adjoint equation is defined as 

= JJq (F + Yt n d Q (59) 

for all kinematically admissible displacements X. As in the previous 
cases (membranes and bending beams), Eq. (59) is identical to Eq. (57) 
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if X = z. Using the adjoint variable method of design sensitivity 
analysis gives 



where 

given 

and 


JJ Q {2yz - Et 2 [ z }} z 22 + 2 Vz 1 1 z 22 + 2(1 " v)z 12^ /4 (1 " ^ } 6t 
2 

If [ 2 yz - l a i ^(z)e i ^(z) ] 6 t dfi (60) 


»«<z> 

and e i ^(z) 

are the stress 

and strain of the extreme fiber, 

as 





Ifii 

2 


(61) 

11 

Et 

Ki + 


0 = ~ 

2(1 - v 2 ) 


22 

Et 

t*22 + «ll) 

(62) 

O — ““ 

2(1 - v 2 ) 

a 12 - - 

Et 



O — 

2(1 + v) 2 

12 



To numerically integrate Eqs. (57) and (60), a one-point Gauss 
quadrature formula on a triangular element is used to correspond to the 
IFAD integration technique for this element. Equations (57) and (60) 
become 

N 

* 7-1 {(F + Yt. )z W}j (63) 

' k=l K 

and 

N 2 

♦7 - l [2yz - l o iJ (z)e 1 :l (z) ]W J 6 t 
k=l i,j-l 


( 64 ) 
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Next consider the functional representing the displacement z at a 

A 

discrete point x as 

i|> g = z(x) = JJ a 6(x - x)z(x)dO (65) 

The first variation of Eq. (65) is 

^8 = Ha 5 ^ x “ x ) z '( x ) dx (66) 

The adjoint equation is defined as 

a u (A,X) = JJ a 6(x - x)A(x)dfl (67) 


for all kinematically admissible virtual displacements A. This equation 

( 8 ) ( 8 ) 

has a unique solution X v , where is the plate displacement due to 

a 

a unit vertical load acting at a point x. Using the adjoint variable 
method of design sensitivity analysis gives 


= ff J y> ( 8) _ pt 2r .(8) .(8) , (8) 

^8 Et L z n X n + z 22 X 22 t 11 X 22 

+ 2(1 - v)z 12 A^ ) ]/4(1 - v 2 )}6t dfl 

= JL [yA (8) - l a ij (z)e ij (A (8) )]6t dfl 
i.J-1 


z A (8) 
Z 22 A 11 


( 68 ) 


The same numerical integration procedure is used as in the 
compliance case. Equation (68) becomes 


*8 


l [yx (8) 


l a ij (z)E ij (A (8) )]w J 6t. 
i,j=l 


(69) 


Finally, consider the functional representing a locally averaged 


stress in the plate as 
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*9 = J7fl «(«(*> V p 


(70) 


where g(a(z)) may be principal stress, von Mises' stress, or some other 

material failure criteria and Up is a characteristic function defined on 

a finite element 0 as 
P 


m 

P 



7^ 


X e n 

p 


X ^ fi 

p 


The first variation of Eq. (70) is 


(71) 


*9 


- // 


Q 


2 

[ l 


3 a' 


( z ' ) ]m d ft 
,1J v J n 


The adjoint equation is defined as [1] 


(72) 


a u (X,X) 


- If 


a 


2 

[ 1 




8 o' 


ij 


a J (X)]m dft 


(73) 


for all kinematically admissible virtual displacement X. Using the 
adjoint variable method of design sensitivity analysis gives 


*9 = //« [ yX<9) " l a ij (z)e i:l (X (9) )]6t d ft + 6t dft 

i,j=l 

(74) 

(9) 

where X is the solution of Eq. (73). For principal stress the last 


terra on the right of Eq. (74) becomes 



//a K’ 11 + ° 22 + 2 w) /<2t) ]" P di! 


(75) 
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where t is defined in Eq. (27). For von Mises' stress, this term 
max 

becomes 


rr 1 U 11x2 11 22 ^ , 22.2 0 , 12.2 xl/2 . ^ 

J/fl- [(a ) - a a + (a ) + 3(a ) ] m p «St da 

(76) 


The IFAD thin shell element is the element that is used for plate 
bending. The membrane effects can be eliminated so that only the 
bending term exists. This element is a hybrid element which uses a 
derivative smoothing technique [4]. The IFAD code evaluates the normal 
and shear stresses at the centroid of the triangle, but Ref. 4 
stipulates that the stresses at the midside nodes of the triangle give 
the most accurate results. The integration of a function <f> from its 
mid-stresses is [4] 


J ♦ dA = | [♦Co, 7 » 7 ) + <H Y » 7 .0) + *(| , \ ,0)] 

(77) 

where A is the area of the triangle and (0, 7 , 7 ), ( 7 , 0 , 7 ) 
and ( j , y ,0) are the area coordinates of the mid-side nodes. To 
achieve the most accurate atresses possible, this numerical integration 
technique is applied on Eqs. (70), (73), and (74). These equations 
become 


- I g(<l n (z))n. f 

n=l 

J7 fl [ I o ij (X)]m dfi 

i, j=l do 2 p 



n n n 


( 78 ) 
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[E] [B] in 4 d - F T d 
n p J 

and 


( 79 ) 


*9 - X < ? K 9> - ! 

k=l n=l i,j=l 

(80) 

where subscript n is the counter for the midside value. The term 
(3g/3t) for principal and von Mises' stress is 


[[a 11 + a 22 + 3r )/(2t, )]m 
L v n n max ' k J p 


(81) 


and 


1 c t 1K2 11 22 ^ , 22 n 2 , 12x2! 

— (a ) " o a + (a ) + 3(a ) ] 

t- L n nn n n J 

k 


1/2 


m 


(82) 


respectively. 

To numerically calculate the adjoint load, the shape function of 
the element must be known so that o^(A) can be calculated using the 
finite element technique described in Section 2.1.1, Eqs. (23) and (25). 


2.2 Calculation Procedure for a Built-Up Structure 
The foundation of the built-up structure design sensitivity 
analysis method is the structural component analysis developed in 
Section 2.1. A built-up structure consists of various structural 
components that interact with each other. This interaction is taken into 
account by generalizing the individual components such that twisting, 
bending, transverse shear terms, etc. are included in the formulation of 
the sensitivity vector. Coordinate system precautions need to be taken 
to insure that the constraints and the sensitivity vectors are 
calculated correctly. In general, if the calculations are performed at 
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the local element coordinate system level, there will be no problem when 
components are oriented differently in the global coordinate system. 

An example of a built-up structure using the structural components 
developed in the previous section would be a bending beam and plate 
problem, where a framework of beams could act as the supporting 
structure for the plate; i.e., a roof structure. 

Figure 6 shows a built-up structure that has design variables 
u=[b(x),h(x),t(x)]^, where b(x) is the width of the beams, h(x) is the 
heigth of the beams, and t(x) is the thickness of the plates. 

It is assumed that the plates are welded along the length of the 
beams. This would infer that the appropriate components in the finite 
element analysis must be chosen to insure kinematic compatibility along 
the component boundaries. 

The energy bilinear form of the system equation is just the sum of 
the plate and beam energy bilinear equations, with an additional beam 
torsion term given as 


A 

a (z,z) = 7 //j)(u)fz z + z z + ftz z + z z ) 
u L JJ sr 1 xx xx yy yy L yy xx xx yy 


+ 2(1 - v)z z ld£l + 7 E, 7-— z z dx a 
xy xy J L J 0 b 12 xx xx 1 



i 


l a u <2 - 2 >plate + 1 a „ <2 ' 2) baam + l k GJ 


xy 


z dx, 
xy Z 


( 83 ) 
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Figure 6. Roof Structure 


where a u ( z > z )pi a te an< * a u^ z,z ^beam are t * ie ener Sy bilinear forms of the 
plate, Eq. (54) and beam, Eq. (31), respectively, G is the modulus of 
rigidity, J is the torsional moment of inertia of the beam, and 
£ represents the local beam coordinate system, where x^ runs along the 
length of the beam. 

In Eq. (83) z represents the beam torsion term. Since each 
xy 

structural component needs to be solved individually to make the 

analysis feasible with an existing finite element code, a relationship 

between z and the beam rotation has to be used. This kinematic 
xy 

compatibility states that for the beam and plate system the term has to 
be equivalent to the relative angle of twist over an element length. 
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That is to say, 


“ q1x 


xy 


(84) 


2 1 
where 0 x^ is the local element rotation at node 2 of the beam, 9 x^ 

is the local element rotation at node 1 of the beam, and L is the beam 

element length* 

The load linear form of the system is 


*U Cz) = //a K + V + Y b bh ^ da (85) 

where Fp is the externally applied plate pressure, y^ and y^ are the 
material densities for the plates and beams respectively, and z is the 
virtual displacement. The state equation is [1] 

a u (z,z) = Z y (z) (86) 

for all kinematically admissible virtual displacements z. 

Since the energy bilinear form of the system equation is just the 
addition of each structural component's energy bilinear forms, the 
design sensitivity equation of the system turns out also to be an 
additive process. The generalized design sensitivity of the built-up 
structure is 

i|»' » ^ 6b + 5h + i|)£ 6t (87) 

The only necessary step to calculating this value is the reformulation 
of the beam to include the torsional term. The energy bilinear form of 
the beam component becomes 

a (z,z) * E(bh 3 /12)z z dx + GJ z 2 dx 

u * J 0 xx xx J 0 xy 


( 88 ) 
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where Eq. (84) defines z X y» 

If the compliance, displacement, and stress functionals of the beam 
- Eqs. (34), (40), and (45), respectively - remain the same, the only 
additional term in the design sensitivity analysis [1] is due to the 
differentiation of the torsional terras in the energy bilinear equation 
with respect to the design variables b and h. This terra is defined as 


/^ [ 6h + 6b]G z X dx 
J 0 L 3h 3b J xy xy 


(89) 


For a beam with a rectangular cross section [5] 


J = bh 3 [ \ - 0.21 (b/h) (1 - )] 

O 1 V 


(90) 


12h 


and the derivatives with respect to the design variables are 


|i-f- 0.042b (h 2 ^ 


(91) 


and 


= bh 2 - 0.42 b 2 (h - ) 

dD 12h 


(92) 


The compliance sensitivity of Eq, (37) becomes 


*4 = % l~ 2yhz ~ (Eh 3 /12)(z xx ) 2 - -f£ G(z xy ) 2 ]dx 6b 
+ [-2ybz - (3Ebh 2 /12)(z xx ) 2 - G(z xy ) 2 ]dx 6h 


( 93 ) 
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the displacement sensitivity of Eq. (43) becomes 


K = /n ["hyX (5) - fr- z * (5) - G z X (5) ]dx 5b 
T 5 J 0 L 12 xx xx 3b xy xy J 


(94) 


+ t [-2YbX <5) - - f S z X <5) ]<!, 5h 

■'ll L 12 xx xx 3h xy xy J 


and the stress sensitivity of Eq. (49) becomes 


I, rL r , v ,(6) Eh .(6) 3J _ ,(6)i. 

M " Jr\ L“nYX -rz— z X - -rr- G z X dx 5b 
6 •'O L 12 xx xx 3b xy xy J 


. fL r_ I F r ™ - Wl< 6 > - 3Ebh „ 1 
■*0 L 2 E z xx p brX 12 Z xx X 


( 6 ) 

xx 


(95) 

$0. x (6 > 

3h xy xy 


]dx5h 


when the torsional term is added, where X^^and are the solutions 

to the adjoint Eqs. (42) and (48), respectively. 

Caution has to be taken when the constraint functionals on the 

system are defined. When a von Mises* stress functional is specified for 

a particular plate element, the allowable beam bending stress term 

- fn t E z m d 6h must be removed from Eq. (95), so that the design 
■'02xxpx ° 

sensitivity of Eq. (87) is calculated only for the von Mises* stress 
functional. In the same way, if an allowable beam bending stress 
functional is specified for a particular beam element, the von Mises' 
stress term 3g/ 3t of Eq. (76) must be removed from Eq. (74), so that 
the design sensitivity of Eq. (87) is calculated only for the allowable 
beam bending stress functional. The application of this procedure is 
shown below for a von Mises* stress functional, a compliance functional 
and a displacement functional on a plate element and an allowable 
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bending stress functional on a beam element of the built-up structure 
of Fig. 6. 

The functional representing a von Mises' stress constraint on a 
plate element is 


.2 1 l/2 


■ //_ [o - a a + <j + 3t ] 

T 10 ■'•'S2 L xx xx yy yy xy J p 


m d Q 


(96) 


The adjoint load linear form is defined as 




i,j=l 3a 

_3g_ 


3a 


JSL 


is_ 


XX 


3a ’ 3 t 
yy xy 


][E] [B] m p dfl d 


(97) 


T - 
F d 


where 3g/3a^ is defined in Eq. (30) and F is the adjoint equivalent 
nodal force. Using Eq. (87) gives the design sensitivity of the roof 
structure as 


+i 0 = I If a [V X 


(10) 


! a 1J (z)c 1 J(X <10) ) + -f 1U. 

i,j“l 


* l f L 0 [-1>TX <10) 

♦ l Jo [-™ (10) 


Eh 3 (10) 
12 Z xx A xx 


3Ebh 2 (10) 
12 Z xx A xx 


G z X (10) ]dx 5b 
3b xy xy J £ 


-|r- G z X (10) ]dx 6h 
3h xy xy J % 

(98) 


where subscript l refers to the beam local coordinate system and 
3g/3t is defined in Eq. (76). 

The functional representing a compliance constraint on a plate 


element is 
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*11 " Hq (F + Yt + Ybh)z dfl (99) 

where yt and ybh are the self weight of the plate and beam, 
respectively. The adjoint load linear form is defined as 

a u (A,X) = J/ fl (F + yt + Ybh) A dfl (100) 

Since Eq. (100) is identical to Eq. (99) if X ■ z, the adjoint equation 
does not need to be solved. Using Eq. (87) gives the design sensitivity 
of the built-up structure as 

♦j, = l JJo[ 2 Yz - l o i ^(z)e i ^(z) ]6t dft 
i,j=l 

♦ 1 fo - e " 3 / 12 <^> 2 - f G < V 2 ]^* «» 

+ I HyKz - (3Ebh 2 /12)(z zx ) 2 - -g G(z x5 ,) 2 ]dx I 6h 

( 101 ) 

where subscript i refers to the beam local coordinate system. 

The functional representing the displacement z at a discrete 
* 

point x is 

i|> 12 = z(x) = Jf a 6(x - x)z(x)d£2 (102) 

The adjoint load linear form is defined as 

a (A, A) = ff n 6(x - x)A(x)dft (103) 

for all kinematically admissible virtual displacements A. This equation 

( 12 ) ( 12 ) 

has a unique solution A , where A is the plate displacement due 
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to a unit vertical load acting at a point x. Using 

Eq. (87) gives the design sensitivity of the built-up structure as 


*; 2 - l //„ [yx < 12) - l 0 i3 (z> £ 13 u (12> )] 
12 “ i,J-l 


+ I t I-hYX <12) - W- z X <12) - £ Gz X (12, ]d* ( 5b 
4 Jq L 12 xx xx 3b xy xy J A 


+ l /o [-W 


(12) _ 3Ebh (12) 

12 z xx xx 


3J „ , (12) i „ 

-rr Gz X dx. 61 
3h xy xy J A 

(104) 


where subscript A refers to the beam local coordinate system. 

The functional representing an allowable bending stress on a beam 
element is 


+13 ■-•<0 2 “*= m p dx 


(105) 


The adjoint load linear form is defined as 


/p y hEX^ m p dx - /q y ^E[B]m p dx d - F T d (106) 

where F is the adjoint equivalent nodal force. Using Eq. (87) gives the 
design sensitivity of the built-up structures as 


*{3 = I JJ n [YX (13) - l a ij (z)e ij (X (13) )]6t dfl 


i.J-1 

3 


, y f L r . v .(13) Eh_ (13) W. . X (13)i. x - 

+ 1 J 0 L-hYX - YT Z xx X xx " lb G Z xy X xy J dx A 615 
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* I ft [-w 


(13) _ 1 _ m 3Ebli .(13) 
2 E2 xx m p 12 Z xx X xx 


tt- G z X (13) ]dx 5h 
3h xy xy J & 


(107) 


where subscript l refers to the beam local coordinate system. 
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CHAPTER III 
NUMERICAL EXAMPLES 

The design sensitivity of a constraint functional ij/ is the 
differential f'of the constraint functional. In order to make sure the 
design sensitivity is accurate, an approximation AiJ> is made using the 
finite difference method. It is very important that an appropriate 
perturbed design variable 5u is selected. If 6u is too small, the 
change in the constraint functional A\|> may be inaccurate due to losses 
of significant digits. If 5u is too large, A\J» will be influenced by 
nonlinearities in the constraint functional, which in turn will cause an 
inaccurate design sensitivity prediction. 

In all the following examples, perturbations in design of 1% and/or 
5% are used, with the exception of the plate bending and built up 
structure examples. Because of the nonlinearity characteristics of 
built-up structural response, the perturbation of 0.1% was used. 

3.1 Membranes 

The finite element membrane model in Fig. 7 is a simple plane 
elastic solid that is restrained at one end and loaded with a 
distributed tensile load at the other end. It contains 80 isoparametric 
elements (IFAD plane stress element type 1104), 289 nodal points, and 
560 degrees-of-freedom, with the design variable being the variable 
thickness u = h(x). The material property constraints. Young's modulus, 
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and Poission's ratio are given as E = 3x10^ psi and v = 0.3, 
respectively. Each finite element is discretized so that a constant 
thickness of h = 0.5 in. is used in order to simplify the sensitivity 
calculation. 

The compliance sensitivity results are shown in Table 1 , where 
A^i * ipj (h + 6h) - i)>(h) and f' and is the predicted value calculated 
from Eq. (9), with design perturbations of 5h * O.Olh and 6h=0.05h. 

The percent accuracy of the sensitivity prediction is calculated using 
* x 100/AiJ^. 

Table 1 . Membrane Design Sensitivity Check for Compliance 


6h 

fyOO 

ipj'Ch+fihj 

♦i 

xl00/ 

O.Olh 

265.302 

262.676 -2.627 

-2.653 

101 .0 

0.05h 

265.302 

252.668 -12.632 

-13.265 

105.0 


Several discrete points shown in Fig. 6 are selected to check 
accuracy of the design sensitivity of the displacement functional of 
Eq. (14). In order to calculate this equation, the strain due to 
the adjoint load is needed. Since the adjoint load is just a unit point 

A 

load at point x, acting in the direction of the displacement, a restart 
of the finite element analysis is all that is needed. For every node 
direction, there is a separate load case that produces a strain 

, which in turn is used to calculate sensitivity of displacement. 
Design sensitivity predictions and differences, with Sh = 0.05h are 
given in Table 2. 
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Table 2. Design Sensitivity Check for Displacement 


Node % 


No. 

Dir 

^(h) 

4» 2 (h+5h) 

^2 

♦*2 

Accuracy 

23 

xl 

2.974E-04 

2.832E-04 

-1.416E-04 

-1.487E-05 

105.0 

27 

xl 

4.058E-04 

3.865E-04 

-1 .932E-05 

-2.029E-05 

105.0 

27 

x2 

-2.248E-04 

-2. 14 IE-04 

1.071E-05 

1.124E-05 

105.0 

185 

xl 

3.298E-03 

3.141E-03 

-1.571E-04 

-1.649E-04 

105.0 

187 

xl 

3.299E-03 

3.142E-03 

-1.571E-04 

-1.650E-04 

105.0 

187 

x2 

-2.014E-04 

-1.918E-04 

0.959E-05 

1.007E-05 

105.0 

365 

xl 

6.633E-03 

6.317E-03 

-3.158E-04 

-3.316E-04 

105.0 

369 

xl 

6.633E-03 

6.317E-03 

-3.158E-04 

-3.316E-04 

105.0 

369 

x2 

-4.000E-04 

-3.809E-04 

1.905E-05 

2.000E-05 

105.0 


To check the stress constraint sensitivity of Eq. (22), the 
equivalent nodal force of the adjoint load on the right of Eq. (25) has 
to be calculated so that e^(X^^) is known for each constrained 
element. This is accomplished by a restart of the original IFAD model, 
with each adjoint load for an element being a separate loading case. 
Design sensitivity results for principal and von Mises' stress 
functionals are given in Table 3 for several finite elements. The 
perturbations are 6h ■ O.Olh and Sh = 0.05h for von Mises' stress and 
6h ■ 0.05h for principal stress. 

The design sensitivity calculation is performed using double 
precision accuracy. Since the finite difference approximation in Tables 
1, 2, and 3 are no smaller than two significant digits of the actual 
constraint functional, loss of significant digits is minimal. 

Therefore, the design perturbation is not to small. With all three 






Table 3. Design Sensitivity Check for Stress 
(a) von Ml8es' Stress with 6h - O.Olh 


El. 

No. 

t 3 (h) 

* 3 (h + 6h) 

4*3 

(^/A* 3 x 100)Z 

1 

9888.882 

9790.973 

-97.910 

-98.889 

101.0 

10 

9989.932 

9891.022 

-98.910 

-99.899 

101.0 

20 

9999.982 

9900.972 

-99.010 

-100.000 

101.0 

21 

8752.850 

8666.188 

-86.662 

-87.528 

101.0 

30 

10024.586 

9925.333 

-99.253 

-100.246 

101.0 

40 

9999.853 

9900.845 

-99.008 

-100.000 

101.0 



(b) von Mises* Stress with 6h 

- 0.05h 


El. 

No. 

* 3 oo 

* 3 <h + 6h) 

A* 3 

* 3 (*3/6*3x100)2 

1 

9888.882 

9417.984 

-470.899 

-494.444 

105.0 

10 

9989.932 

9514.222 

-475.711 

-499.497 

105.0 

20 

9999.982 

9523.793 

-476.189 

-499.998 

105.0 

21 

8752.850 

8336.048 

-416.802 

-437.642 

105.0 

30 

10024.586 

9547.225 

-477.361 

-501.230 

105.0 

40 

9999.853 

9523.670 

-476.183 

-499.993 

105.0 



(c) Principal Stress with 6h 

- 0.05h 


El. 

No. 

* 3 <t») 

* 3 (h + 6h) 

4*3 

*3 

(*^/A* 3 x100)Z 

1 

10582.770 

10078.829 

-503.941 

-529.138 

105.0 

10 

9987.061 

9511.487 

-475.574 

-499.353 

105.0 

20 

10000.013 

9523.823 

-476.191 

-500.001 

105.0 

21 

9660.279 

9200.266 

-460.013 

-483.014 

105.0 

30 

10012.976 

9536.617 

-476.808 

-500.649 

105.0 

40 

9999.987 

9523.797 

-476.189 

-500.000 

105.0 
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constraint functionals, the design sensitivity results compared to the 
finite difference approximation are excellent. This infers that the 
design perturbation is not too large as to cause significant 
nonlinearity effects in the calculation of design sensitivity. 

It is interesting to note that in Tables 1, 2, and 3 the finite 
difference approximation is nearly 1% of the constraint functional 
when 6h = O.Olh and nearly 5% of the constraint functional when 
h = 0.05h. The results also show that as 6h approaches zero, 
approaches one. 


3.2 Bending of Beams 

A cantilevered beam finite elment model shown in Fig. 8 is loaded 
with a constant distributed force f(x) = 0.03 Ib/in. along the entire 
length of the beam. It contains 20 IFAD beam elements of type 0501, each 
3" in length, 121 nodal points, and 40 degrees-of-f reedom, with design 
variables u = [b(x),h(x)] , the width and height of the beam. In 
Fig. 8, the element numbers are along the top of the beam and the node 
numbers are along the bottom of the beam. The beam has a rectangular 
cross-section with constant width and heigth b = 0.5 in. and h = 0.75 
in., respectively. This gives the moment of inertia Iy = 0.01758 in.^ 

O 

and the cross-sectional area A c = 0*375 in . The material property 

7 

constants, Young’s modulus and Poisson’s ratio are E = 3x10 psi and 
v = 0.3, respectively. Self weight is included in the analysis. 
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Figure 8. Cantilever Beam Finite Element Model 

The compliance sensitivity results are shown in Table 4, where 
Ai|» 4 = i|» 4 (u + 6u) - \j> 4 (u) and ip 4 is the predicted value calculated from 
Eq. (39), with design perturbations of 6b = 0.05b, 6h = 0.05h and 6b = 
0.01b, 6h = O.Olh. 

Table 4. Beam Design Sensitivity Check for Compliance 


6h 

6b 

V u) 

ij» 4 (u+6u) 

^4 

*4 

% 

Accuracy 

0.05h 

0.05b 

1.3684 

1.3153 

-0.0531 

-0.0601 

113.1 

O.Olh 

0.01b 

1.3684 

1.3566 

-0.0118 

-0.0120 

102.0 


These results show that nonlinearities in the compliance of the 
cantilevered beam are highly evident. If too large a perturbation in 
design is chosen, the sensitivity will be inaccurate. 
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Several discrete points along the beam are selected to check the 
accuracy of design sensitivity of the displacement functional of 
Eq. (44). In order to calculate this equation, the beam curvature due 
to the adjoint load is needed. Since the adjoint load is just a unit 

A 

point load at the point x acting in the -z direction, a restart of the 
finite element analysis is all that is needed. Displacement results are 
shown in Table 5 for design perturbations of 5b « 0.05b, 6h = 0.05h 
and 6b * 0.01b, 6h = O.Olh. 

In both Tables 5(a) and 5(b), results show that the design 
sensitivity compared to the finite difference approximation is good, 
with the exception of node 3. Since the finite diffenrence 
approximation is not too small in comparison to the constraint 
functional, loss of significant digits is not a valid reason for this 
inconsistency. The accuracy decreases as the node location approaches 
the restrained end of the beam. This is most likely due to the 
restraining effect, which causes rigidity in the beam and in turn gives 
smaller deflections. 

To check the stress constraint sensitivity of Eq. (52), the 
equivalent nodal force of the adjoint load on the right side of Eq. (51) 
has to be calculated so that the curvature X^^ is known for each 
constraint element. This is accomplished by a restart of the orginal 
IFAD model, with each adjoint load for each element being a separate 
loading case. Allowable bending stress results for several finite 
elements are shown in Table 6 for a design perturbation of 5h = 0.05h 


and 6b = 0.05b. 



Table 5. Beam Design Sensitivity Check for Displacement 
(a) 6b= 0.01b and 6h = O.Olh 
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Node 


ij> 5 (u+6u) 


*5 


< 

o 

o 

X 

3 

7.8353E-03 

7.6473E-03 

-1.8804E-04 

-1.6082E-04 


85.5 

6 

4.4162E-02 

4.3102E-00 

-1.0604E-03 

-9.9992E-04 


94.2 

9 

1.0181E-01 

9.9365E-02 

-2.4451E-03 

-2.3554E-03 


96.3 

12 

1.7317E-01 

1.6901E-01 

-4 . 1590E-03 

-4.0441E-03 


97.2 

15 

2.5229E-01 

2.4623E-01 

-6.0594E-03 

-539211E-03 



18 

3.3493E-01 

3.2686E-01 

-8.0447E-03 

-7.8836E-03 


97.9 

21 

4.1857E-01 

4.0852E-01 

-1 .0054E-02 

-9.8701E-03 


98.2 




(b) 6b ■ 0.05b and 6h = 0 

.05h 


Node 

’I'jCu) 

t 5 (u+6u) 


*5 

i(^xl 00 AiJ> 5 

3 

7.8353E-03 

6.9743E-03 

-8.6100E-04 

-8.0412E-04 

93.4 

6 

4.4162E-02 

3.9306E-02 

-4.8555E-03 

-4.9959E-04 

102.9 

9 

1.0181E-01 

9.0619E-02 

-1.1120E-02 

-1 . 1777E-02 

105.2 

12 

1.7317E-01 

1.5412E-01 

-1.9043E-02 

-2.0220E-02 

106.2 

15 

2.5229E-01 

2.2454E-01 

-2.7745E-02 

-2.9605E-02 

106.7 

18 

3.3493E-01 

2.9810E-01 

-3.6835E-02 

-3.9418E-02 

107.0 

21 

4.1857E-01 

3.7253E-01 

-4.6034E-02 

-4.9350E-02 

107.2 



Table 

6. Beam Design 

Sensitivity Check for Stress 

El. 

V u) 

i|» 6 (u+6u) 

% 

4*5 ^x100/Ai|> 6 

1 

4932.767 

4609.263 

•323.504 

-349.234 108.0 

5 

3110.255 

2906.201 

■204.055 

-219.226 107.4 

10 

1420.622 

1327.375 

• 93.287 

- 99.103 106.2 

15 

385.008 

359.664 

25.344 

- 26.076 102.9 

20 

3.295 

3.067 

• 0.228 

- 0.147 64.7 
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Near the end of the beam, where the allowable bending stress is 
near zero, the design sensitivity decreases significantly. Since the 
design sensitivity is the derivative of the constraint functional, and 
the derivative is physically interpreted as the slope of the beam, this 
decrease can be attributed to the large increase in slope at the free 
end of the beam. 


3.3 Bending of Plates 

The clamped plate finite element model shown in Fig. 9 is uniformly 
loaded with a pressure f(x) - -1.5 lb/in. in the z direction. Since the 
model is symmetric along two planes, only one quarter of it needs to be 
analyzed and symmetric boundary conditions need to be applied. The 
quarter model contains 100 IFAD triangular thin shell elements of type 
1601, with only the bending terms active. It has 61 nodal points and 
140 degrees-of-freedom. 

The design variable is the plate thickness u = t(x), and the 
material property constants. Young’s Modulus and Poisson’s ratio, are 
E 30.5x10^ psi and v = 0.3, respectively. The constant plate thickness 
is t = 0.4 in. and the self-weight of the plate is neglected. 

Compliance sensitivity results are shown in Table 7, where 

= \|)^(t + 6t) - ^y(t) and is the predicted value calculated from 
Eq. (60), with design perturbations of 6t = O.Olt and <5t = 0.05t. 

Both perturbations for the compliance constraint functional give 
good correlation between the design sensitivity and the finite 
difference approximation. This implies that a five percent change in 
thickness is acceptable when making design improvements. 
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Xg • 


Figure 9. Bending Plate Finite Element Model 
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Table 7. Plate Design Sensitivity Check for Compliance 


<5t 

t 7 (t) 

y t+6t) 

Ai |> 7 

*7 

i|> 7 xl00 / Ai|> 7 

.05t 

4.9625 

4.2868 

-0.6757 

-0.7204 

106.6 

.Olt 

4.9625 

4.8165 

-0.1459 

-0.1441 

98.7 


Several discrete points in Fig. 8 are selected to check accuracy of 
the design sensitivity of the displacement functional in 
Eq. (69). In order to calculate this equation, just as in the membrane 
case, the strain due to the adjoint load is needed. A restart 

of the finite element analysis, using a unit point load in the -z 
direction for each selected point as a separate load case, accomplishes 
this task. Some displacement results are shown in Table 8 for a design 
perturbation of 6t = O.Olt. 


Table 8. Plate Design Sensitivity Check for Displacement 


Node 

yo 

y t+6t) 

% 

OO - 

V 8 

xlOO/Aik, 

O 

14 

1.9033E-03 

1.8473E-03 

-5.5976E-05 

-5.2596E-05 


94.0 

15 

3.0497E-03 

2.9600E-03 

-8.9692E-05 

-8.5144E-05 


94.9 

24 

1.9033E-03 

1.8473E-03 

-5.5976E-05 

-5.2596E-05 


94.0 

27 

1.1258E-02 

1.0927E-02 

-3.3109E-04 

-3.2620E-04 


98.5 

31 

9.7772E-03 

9.4897E-03 

-2.8754E-04 

-2.8309E-04 


98.4 

35 

3.0497E-03 

2.9600E-03 

-8.9692E-05 

-8.5144E-05 


94.9 

38 

1.8450E-02 

1.7908E-02 

-5.4262E-04 

-5.3957E-04 


99.4 

47 

1.1258E-02 

1.0927E-02 

-3.3109E-04 

-3.2620E-04 


98.5 

48 

1.8450E-02 

1.7908E-02 

-5.4262E-04 

-5.3957E-04 


99.4 

57 

4.2639E-03 

3.9443E-03 

-1.1952E-04 

-1.1410E-04 


95.5 

60 

2.5084E-02 

2.4366E-02 

-7.3773E-04 

-7.3688E-04 


99.9 

61 

2.6942E-02 

2.6150E-02 

-7.9237E-04 

-7.9223E-04 


100.0 
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The design sensitivity results of Table 8 seem to substantiate the 
fact that nodes near a clamped edge give a somewhat less accurate 
prediction of the design sensitivity. However, these results are still 
considered good, since there is less than 15% deviation from the 
approximate finite difference result. Since the finite difference 
result is only an approximation, it is reasonable to say that the design 
sensitivity prediction is good every where along the plate. 

A number of elements are selected to check the design sensitivity 
for von Mises 1 stress in Eq. (80). Before this an be evaluated, the 
nodal force of the adjoint load on the right side of Eq. (79) has to be 
calculated so that the strain e^(X^^) is known for each constraint 
element. This is accomplished by a restart of the original IFAD model 
with each adjoint load for each element being a separate loading case. 
Von Mises 1 stress results are shown in Table 9 for a design perturbation 
of 6t « O.OOlt. 

Recall that the design sensitivity and the constraint functional 
was calculated using an integration technique where the stresses were 
calculated at the midside nodes of the triangular element. This 
integration technique was used so as to get the best design sensitivity 
results as possible for the finite elements. This is the technique used 
because the element is a hybrid element [4]. The results in 
Table 9 are excellent when this technique is used. These results 
further substantiate the fact that as design perturbation approaches 
zero, approaches one. 
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Table 9. Plate Design Sensitivity Check for von Mises' Stress 


El. No. 'l'g(t) 

l|»g(t+6t) 


*9 

*9 

X 

>—* 

o 

o 

> 

1 

172.1862 

171.8423 

-0.3439 

-0.3445 


100.17 

2 

752.7383 

751.2351 

-1.5032 

-1.5049 


100.11 

3 

172.1861 

171.8423 

-3.4386 

-3.4446 


100.17 

4 

752.7383 

751.2351 

-1.5032 

-5.0495 


100.11 

10 

1442.8840 

1140.0025 

-2.8815 

-2.8841 


100.09 

17 

2457.3522 

2452.4447 

-4.9074 

-4.9125 


100.10 

25 

1149.0097 

1146.7151 

-2.2947 

-2.2970 


100.11 

26 

1331.9565 

1329.2966 

-2.6600 

-2.6624 


100.10 

27 

1149.0097 

1146.7151 

-2.2946 

-2.2970 


100.11 

28 

1331.9565 

1329.2966 

-2.6600 

-2.6662 


100.09 

34 

973.8975 

971.9526 

-1.9450 

-1.9465 


100.08 

44 

1442.8840 

1440.0025 

-2.8815 

-2.8841 


100.09 

49 

1284.1558 

1281.5913 

-2.5645 

-2.5667 


100.08 

50 

1278.5066 

1275.9534 

-2.5532 

-2.5554 


100.09 

51 

1274.1558 

1281.5913 

-2.5645 

-2.5667 


100.08 

52 

1278.5066 

1275.9531 

-2.5532 

-2.5554 


100.09 

57 

1110.0286 

1107.8118 

-2.2168 

-2.2188 


100.09 

68 

973.8975 

971.9526 

-1.9449 

-1.9465 


100.08 

73 

1470.1247 

1417.2887 

-2.8360 

-2.8384 


100.08 

74 

1629.1926 

1625.9391 

-3.2536 

-3.2563 


100.08 

75 

1420.1247 

1417.2887 

-2.8360 

-2.8384 


100.08 

76 

1629.1926 

1625.9391 

-3.2536 

-3.2561 


100.08 

83 

2457.3522 

2452.4447 

-4.9074 

-4.9125 


100.10 

91 

1110.0286 

1107.8118 

-2.2168 

-2.2188 


100.09 

97 

1884.4503 

1880.6870 

-3.7633 

-3.7660 


100.07 

98 

2010.1163 

2006.1021 

-4.0143 

-4.0170 


100.07 

99 

1884.4503 

1880.6870 

-3.7633 

-3.7660 


100.07 

100 

2010.1663 

2006.1021 

-4.0143 

-4.0170 


100.07 
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3.4 Built-Up Structure 

A built-up structure that uses both beams and plates is shown in 
Fig. 10. Since it is symmetric along two planes, only a quarter is 
modeled. The built-up structure is clamped on two edges, with symmetric 
boundary conditions applied along the other two edges. The plates and 
the beams are considered to be welded. A uniform pressure f(x) = -1.5 
lb/ in. is applied on the top surface of the plates. The model contains 
100 IFAD triangular thin shell elements of type 1601, with only the 
bending terms active, and 20 IFAD beam elements of type 0501. There are 
61 nodal points and 140 degrees-of-freedom. There are three design 
variables, beam width, beam height, and plate thickness. The material 
constants, Young's modulus and Poisson's ratio, for both the beams and 
the plates are E = 30.5x10^ psi and v = 0.3, respectively. Self weight 
is neglected and the initial design variables are b = 0.5 in., 
h = 0.75 in., and t ■ 0.4 in. 

Compliance sensitivity results are shown in Table 10, where 
A^n = ^(u + “ 'I'jqCu) anc * ^11 * s t * ie P re dicted design sensitivity, 

using Eq. (101) with a 1% design perturbation for all the design 
variables. 

Several discrete points in Fig. 10 are selected to check the 
accuracy of design sensitivity of the displacement functional in 
Eq. (104). To calculate this equation, the strain e^(X)^^ due to the 
adjoint load is obtained by doing a restart of the finite element 
analysis. A unit point load in the -z direction for each selected point 



Figure 10. Built-up Structure Finite Element Model 
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Table 10. Built-up Structure Design Sensitivity 
Check for Compliance with 6h = O.Olh 
(Sb = 0.01b, and 6t = O.Olt 


ij^Cu+du) ♦Jj \|»'i xlOO/ Ai|>j j 


3.7200 3.6017 -0.1183 -0.1189 100.5 


is specified as a separate load case. The design sensitivity results 
for these selected points are shown in Table 11, where a 1% perturbation 
for each design variable is used. 


Table 11. Built-up Structure Design Sensitivity 

Check for Displacement with 6h = O.Olh, 
6b = 0.01b, and 6t = O.Olt 


Node 

♦j^u) 

^12^ U+ &0 

^12 

CM 

rH 

^lOO/Ai^ 2 

13 

4.768E-04 

4.616E-04 

-1.520E-05 

-1.453E-05 

95.6 

15 

2.294E-06 

2.222E-03 

-7.283E-05 

-6.750E-05 

92.7 

17 

3.064E-03 

2.966E-03 

-9.723E-05 

-9.152E-05 

94.1 

25 

4 . 125E-03 

3.994E-03 

-1.313E-04 

-1.255E-04 

95.6 

27 

5.421E-03 

8 . 153E-03 

-2.680E-04 

-2.621E-04 

97.8 

35 

2.294E-03 

2.221E-03 

-7.283E-05 

-6.750E-05 

92.7 

37 

1.095E-02 

1 .060E-02 

-3.484E-04 

-3.448E-04 

99.0 

39 

1.485E-02 

1.438E-02 

-4.724E-04 

-4.737E-04 

100.3 

47 

8.421E-03 

8.153E-03 

-2.680E-04 

-2.621E-04 

97.8 

49 

1.755E-02 

1.699E-02 

-5.580E-04 

-5.621E-04 

100.7 

55 

1.955E-02 

1.892E-02 

-6.215E-04 

-6.296E-04 


57 

3.064E-03 

2.966E-03 

-9.723E-05 

-9 . 152E-05 

94.1 

59 

1.485E-02 

1.438E-02 

-4.724E-04 

-4.737E-04 


61 

2.025E-02 

1.961E-02 

-6.439E-04 

-6.535E-04 

101.5 


Both the compliance and displacement design sensitivity results 
show good correlation with the approximated finite difference. This 
indicates that the method of design sensitivity is a good method for 





53 


predicting the response of a built-up structure for these constraint 
functionals. This is substantiated by the fact that correlation between 
the design sensitivities and the finite difference approximations for 
the built-up structure in Tables 10 and 11 are not that different than 
those in Tables 7 and 8 where plate bending results are shown. 

To check the stress constraint sensitivity of Eq. (98), the 
equivalent nodal force of the adjoint load on the right of Eq. (97) has 
to be calculated so that e^(X)^^ is known for each constrained 
element. A restart of the original IFAD built-up structure model is 
made, with each adjoint load for a finite element being a separate 
loading case. The design sensitivity results for the von Mises' stress 
functional are given in Table 12, with design perturbation 6t = O.OOlt. 

The design sensitivity results compared to the finite difference 
approximations fluctuate more than would be expected, considering the 
excellent correlation of the von Mises' stress design sensitivites for 
the bending plate in Table 9. It doesn't appear that the problem is 
caused by the design perturbation being too small, because the finite 
difference approximation is no smaller than three significant digits of 
the actual constraint functional. It is possible that the nonlinear 
response of the built-up structure is partially the cause of the 
inconsistencies in the design sensitivities, but not the whole problem, 
since the design perturbation is quite small. 

It is most likely that the beam contribution to the built-up 
structure is the major influence. The beam seems to be more influenced 
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by nonlinearities, as shown in Tables 4, 5, and 6. It is most probable 
that if a more complex beam element, such as a cubic beam, was used, 
design sensitivity results would improve. Unfortunately the IFAD code 
does not currently support this type of beam, so this cannot be verified 
in this study. 


Table 12. Built-up Structure Design Sensitivity Check 
for von Mises' Stress with 6t = O.OOlt, 

6b = 0.001b, and dh = O.OOlh 


Element 


ij» 10 (t+6t) 

1 

127.8318 

127.5490 

2 

549.7750 

548.5339 

3 

127.8318 

127.5490 

4 

549.7750 

548.5339 

10 

1074.7414 

1072.3202 

17 

1828.5601 

1824.4411 

25 

875.4135 

873.4610 

26 

985.6012 

983.3797 

27 

875.4135 

873.4610 

28 

985.6012 

983.3747 

34 

735.6485 

734.0066 

44 

1074.7414 

1072.3202 

50 

962.8944 

960.7381 

51 

944.5642 

942.4288 

57 

804.5913 

802.7656 

68 

735.6485 

734.0066 

73 

1074.1010 

1071.6982 

74 

1243.6633 

1240.8853 

75 

1074.1010 

1071.6982 

76 

1243.6633 

1240.8853 

83 

1828.5601 

1824.4411 

91 

804.5913 

802.7656 

97 

1401.1703 

1398.0144 

98 

1517.4522 

1514.0511 

99 

1401.1703 

1398.0144 

100 

1517.4522 

1514.0511 


4 +10 

+10 

tj 0 *100/ 

-0.2828 

-0.3339 

118.1 

-1.2411 

-1.2230 

98.5 

-0.2828 

-0.3339 

118.1 

-1.2411 

-1.2230 

98.5 

-2.4212 

-2.7653 

114.2 

-4.1190 

-4.2965 

104.3 

-1.9525 

-2.3201 

118.1 

-2.2215 

-2.6663 

120.0 

-1.9525 

-2.3201 

118.1 

-2.2215 

-2.6663 

120.0 

-1.6419 

-2.2807 

138.9 

-2.4212 

-2.7653 

114.2 

-2.1563 

-2.4538 

113.8 

-2.1354 

-2.7966 

131.0 

-1.8257 

-2.2406 

122.7 

-1.6420 

-2.2807 

138.9 

-2.4028 

-2.7000 

112.6 

-2.7780 

-3.2605 

117.4 

-2.4028 

-2.7060 

112.6 

-2.7780 

-3.2605 

117.4 

-4.1190 

-4.2965 

104.3 

-1.8257 

-2.2406 

122.7 

-3.1559 

-3.2171 

101.9 

-3.4011 

-3.6964 

108.7 

-3.1559 

-3.2171 

101.9 

-3.4011 

-3.6964 

108.7 
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CHAPTER IV 
CONCLUSIONS 


The results of this study indicate that it is feasible to implement 
the theoretical design sensitivity analysis of Ref. 1 with an existing 
finite element code. Calculations of the design sensitivites can be 
accomplished outside of the finite element code using only the 
postprocessing data. In addition, the results show that accurate design 
sensitivity can be predicted without the uncertainty of numerical 
accuracy associated with the selection of a finite difference 
perturbation. However, results also indicate that the integration 
technique used in the calculation of the design sensitivity and the 
knowledge of the exact finite element shape functions used for each 
finite element in the finite element code is important in getting the 
most accurate design sensitivities possible. 

Results of the built-up structure indicate that design 
sensitivities of the built-up structure cannot be any more accurate than 
the accuracy of the individual components. Care must be taken to use 
the best integration techniques and the same shape functions as the 
finite element analysis for the individual components, if at all 


possible 
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PROGRAM SENSIT 

C ft*********************************************** ****** ******* 
CK* 

CP* SENSITJ THE MAIN PROGRAM bOU CALCULAI 1NG L'ESIUN SFNSU 1VITY 


CP* 

CP* ************************************************** ********** 
CP# 

cp* description: 

CP* 

CP* ' SENSIT ' IS THE MAIN PROGRAM TOR THE DESIGN SENSITIVITY 

CP* VECTOR CALCULATION* IT CHECKS THE ACCURACY UP THE 

CP* SENSITIVITY 10 THE FINITE DIFFERENCE METHOD IF THE 

CP* TWO FINITE ELEMtNf ANALYSES ( THE ORIGINAL ANALYSIS 

CP* AND THE PERTURBED ANALYSIS* ) ARE ALREADY CREATED* 

CP* MAX. OF 500 ELEMENTS* IF MORE ELEMENTS ARE NEEDED THE 

CF* SVECTR • MON COMMON BLOCK EILE NEEDS TO CHANGED. 

CP* 

CP************************************************************* 
INCLUDE ' C AEGSDR * INC J 1MPLIC.SPC' 

INCLUDE ' C AEGSDR. INC J CNTL ♦ MON ' 

INCLUDE ' CAEGSDR . INC J SVECTR. MON' 

C 

EQUIVALENCE <NDAT<46> r NELM) r ( ND A T ( 3 B j ) r IPNAME ) 

C 

DIMENSION PS1 B ( 2 ) , IPNAM£(2) r IPNAM1 (2) * IPNAM2C2) 

CHARACTER YESN0*1 
C 
t: 


c 

c 

c 


c 

c 

c 


c 

20 


c 

c 

c 


NT = 0 
ISAC = 0 
LCS ■ 1 
NLC * 1 


ASK FOR FINITE ELEMENT MODEL FILE NAME 
PRINT *f' ' 

PRINT * t ' ENTER THE ORIGINAL PROBLEM NAME <1-8 CHARS)' 
READ < 5? 1009 ) IPNAM1 

ASK FOR CONSTRAINT TYPE 


PRINT * 9 ' ' 

PRINT * * ' ENTER CONSTRAIN! TYPE.* 
PRINT ft' 

PRINT *t' 

READ ( 5 ? 1006 ) 1CT 

IF ( I CT . NE . 3 ) GO TO 20 

PRINT *.' ' 

PRIM! * r ' EN I ER STRESS IYPEJ 1 = 
PRINT *•' 2 = 

PRINT *?' 3 = 

READ < Sr 1006) 1ST 


3 « COMPLIANCE ' 

2 = D1SPL ACfc MEN T ' 

3 = STRESS' 


PRINCIPAL ' 

VON MIS'ES ' 

BEAM ALLOWABLE' 


PRINT *» 'CALCULATING ADJOINT LOADS <Y/N> Y ' 
READ (Zf 1000) YESNO 
I F < YESNO • EW • ' Y ' ) ISAC = 3 
IF( ISAC. EQ. 1 ) GO TO 40 


ASK FOR THE PERTURBED UNITE ELEMEN T NAME 


PRINT *f 
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PRINT t ? 4 ENTER THE PERTURBED PROBLEM NAME <J-8 CHARS)' 

READ ( 5> 1009 ) JPNAM2 
C 

C OPEN AN OUTPUT KILE 
C 

OPEN ( UN I T= 1 0 1 NAME- ' RESULT . DAT ' * T YPE= ' NEW ' ) 

WRITE ( 10 f 1007 ) ICT 
IF(ICT.NE*3) GO TO 40 
WRITE ( 10 * 1008 ) 1ST 
C 

C GET LOAD CASE START AND END NOS • 

C 

40 PRINT ' 

PRINI # t ' EN t ER NO* UF LOAD CASES TO BE PROCESSED 4 
READ < 5 1 1006 ) NLC 

PRINT tf 'ENTER FIRST LOAD CASE NO* 4 
READ (5*1 006 ) LCS 
■50 LCS = LCS - I 

C 

C GET CONSTRAINED ELEMENTS h OK STRESS CALCULATION 
C 

IF(ICT*NE.3> GO TO 90 
PRINI *»' ' 

IFdSAC.EQ. 1) GO TO 60 
PRINT 2000 * NLC 
GO TO 70 

C 

C GET ELEMENIS FOR ADJOINT LOAD CALCULATION 
C 

60 PRINT 2001* NLC 

70 DO SO NE=1 * NLC 

80 READ <5*1 006 ) ICE(NE) 

90 DO 800 NC=1 * NLC 

IF ( ISAC » EG ♦ J ) GO TO 100 
WRITE < 10 * 1010) ICE(NC) 

C 

C GET SENSITIVITY VECTOR 
C 

100 I F < NT « EG ♦ 1 ) GO TO 110 
IPNAME < 1 ) = I PHAM 1 < 1 ) 

IPNAME ( 2 ) = 1PNAMK2) 

GO TO 120 

110 IPNAME ( 1 ) = I PNAM2 ( 1 ) 

IPNAME ( 2 ) « I PNAM2 < 2 ) 

120 CALL GETSEN <PSI B* N f * NELM * IPNAME ) 

IF ( ISAC * EG* 1 ) GO TO 910 
IF < NT . GT * 1) GO TO 200 
GO TO 100 
c: 

C CALC* CHANGE IN PLATE THICKNESS* CHANGE IN BEAM WIDTH AND DEPTH 
C 

200 N I = 0 

DT - DABS< TM<2)-IM(1>) 

DB = DABS<BW< 2 ) -BW < 1 ) ) 

DH * DABS(BH(2)~BH(1) ) 

TD = DABS ( PB < 2 ) -PB < 3 ) ) 

C 

C CALCULATE THE PERCENT ACCURACY 
C 

DPSJ.BDB = 0*0 
DO 300 1=1 r NELM 
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DPSIBDB = DPSIBDB+DPSITU>*PT+KiPSlB<I>*DB 

* 4* DPS I H ( I ) *DHTDPSITP ( ) ) *T H 
300 CONTINUE 

PNUM = DPS 1 PUP* ICO 
400 PDEN * PS IP ( 2 ) ~PS IP ( 1 ) 

IF ( PDEN ♦ EG • 0 > GO TO 801 
PACCUR * PNUM /PDEN 
C 

WRITEdOr*) 7 ' 

WRITE ( 10 f 1004 ) DT t DB t DH » *1 B 
URITEdOv 1005) PDEN 
WRITEdO>1002) DPSIPUP 
WRITE <10? 100J ) PACCUR 

800 CONTINUE 
C 

GO TO 900 
C 

801 PRINT *? 'DPSKPmip =-* 0' 

C 

900 CLOSE (10) 

910 CONTINUE 

C 

C 

1000 FORMAT < A > 

1001 FORMAT < F8 « 5) 

1002 FORMAT dXr ' DPSI ( B ) #DELTAB= ' ? E16 . 8 ) 

1003 FORMATdX, 'PERCENT ACCURACY- ' * F16 . 8 ) 

1.004 FORMAT < 1 X • ' CHANGE IN MEMBRANE IHICKNESS = ' i F8 . f.,/ , IX * 'CHANGE 

* IN BEAM WIDTH = ' r F8 * 5 > / tl X * ' CHANGE IN PFAM DEPTH =',F8.5r 
*/ ? IX * ' CHANGE IN BENDING PLATE THICKNESS *'»F8.5) 

1005 FORMAT dX > 'PSI ( P+DP> - PSI(H) * ',F16.8) 

1006 FORMAT (14) 

1007 FORMAT ( 1X» ' ***C0NSTRA1N f 1 YPF ^'rl4> 

1008 FORMAT (/r IX r ' TRESS TYPE «'r IA> 

1009 FORMAT < 2A4 ) 

1010 F0RMAT(1Xf 'CONSTRAINT ELEMENT IS 'rI4) 

2000 FORMAT ( IX » 'ENTER' flAr ' CONSTRAINT El EMEN fS r FOLLOW EACH BY 

* A RETURN') 

2001 FORMAT ( lXf 'ENTER' > 14 r ' ELEMENTS THAT ARE fU HAUE AN ADJOINT 

* LOAD CALCULATED. FOLLOW EACH BY A RETURN.') 

C 

C 

END 
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SUEiROUl I HE AL 16 < X f Y» C j AL » TU > I HK » 1 T ) 
CF**************************** a******************************** 

CP* 

CP* AL 16: ADJOIN f LOAD CALCULATION FOR TRIANGULAR El k MEN f 1601 
CP* 

c;p***************** ************************ ******************** 

CP* 

i:p* description: 

CF-‘* 

CP* ' AL.1 6' CALCULATES THE ADJOINT LOADS’ FOR* THF 7 R) ANGULAR 

CP* PLATE LENDING ELEMENT* AL = LCJ*LP:i*MP WHERE 

CP* CC3 IS THE DERIVATIVE OF THE STRESS FUNCTION VECTOR 

CP* TIMES THE ELASTICITY MATRIX* SINCE IE AD CALCULATES 

CP* STRESS RESULTANTS FIRST r L AL 3 MUST HE MULTIPLIED PY 

CP* PLATE THICKNESS DIVIDE PY 2* 

CP* CPI IS A 3X9 MATRIX IN COLUMNS 4fS * 6 UF CW3. MP 

CP* IS THE CHARACTERISTIC FUNCTION THAT IS 1/AREA OF THE' 

CP* ELEMENT THAT IS CUNSTKAlNkfD AND ZERO FOR ALL THE OTHER 

CP* ELEMENTS. 

CP* 

CP ********************************************** *************** 

CP* 

CP* X THE LOCAL ELEMENT X COORDINATE 

CP* Y THE LOCAL ELEMENT Y COORDINATE 

CP* C MATRIX CC3 * CDG3* CE 3*T/2 3X3 MATRIX 

CP* AL ADJOINT LOAD VECTOR 

CP* TP TRANSFORMATION MATRIX 

CP* THK MATERIAL THICKNESS 

CP* IT M IDS IDE NODE COLUMN LOCATOR 

CP* 

(:;p******* ****************************************** ************ 

c 

INCLUDE ' L AEGSDR . INC 3 1MPL1C.SPC' 

INCLUDE 'CAEGSDR.INCJ CN7L.MUN' 

C 

EQUIVALENCE < NDA1 ( 1 4 > * i PR ) 

C 

DIMENSION X<3> fY<3> fGPTS<3»3> »XL<6> »YL<6> f W(18f7> fC<3> » 

* F<9> fAL<18> fK(6f3> fTD<6f6> fRP<6f3> 

C 

DATA GPTS/O.ODOf .3D0f *5D0f .SDOfO.ODOf .DUOf .SDOfO.ODO/ 

C 

C*** INITIALIZE VARIABLES 
C 

DO 10 1 = 1 f V 

io f ( i ) = o.ono 

DO 12 1=1 f 18 
12 AL ( 18 ) = O.ODO 

DO 14 1 = 1 f 6 
DO 14 J= 1 f 3 

R<IfJ) = O.ODO 
14 CONTINUE 

C 

AREA = EUTRIA(XfY) 

XNP = 1. DO/ AREA 
C 

e*** GET LOCAL X AND Y COORDINATES 
C 

CALL M0VESP(XLfXf3*IPK> 

CALL MOVESP <YLfY f3#IPR ) 

CALL SF1501 (XLf YLfGF’TS ( 1 f IT) fWfEM) 
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DO 50 M- 1 f V 
GASH = 0 • ODO 
DO AO J“1 * 3 

GASH ■ CASH+L< J>tW(Mf J+3> 

40 CONTINUE 

F ( M ) = F < M ) +GASH*X HP 
50 CONTINUE 

[;** ROTATE ADJOINT LOAD VECTOR I U GL.OPAL COORDINATE SYSTEM 
C 

N = 1 

DO 60 MM=1*3 
R < 3 * MM ) = F < N ) 

R(4?MM) = F <N+1 ) 

R < 5 * MM ) = F ( N+2 ) 

N = N+3 

60 CONTINUE 

CALL UMXAD ClD>RfRP>6>3f6) 

N1 “0 
DO 80 M=l>3 
DO 70 Ml = 1 » 6 

AL < Ml + N1 ) « RP(MlfM) 

70 CONTINUE 

N1 = N 1+6 
80 CONTINUE 

DO 90 N2 - I * 1 8 

AL(N2) =* AL<N2>* IHK/2.0D0 
VO CONTINUE 

C 
C 

RETURN 

END 



FUNCTION AREAU<X>Y> 

CP******************************************************** ***** 

CP# 

cp* areag: calculates the area of a stray oh i « ided four or 

CP* EIGHT NODE ELEMENT. 

CP* 

CP************************************************************* 

CP* 

CP* X GLOBAL X COORDINATES 

CP* Y GLOBAL Y COORDINATES 


CP* 

CP************************************************************* 

C 

INCLUDE ' CAEGSDR. INCD IMPLIC.SPC' 

INCLUDE ' CAEGSDR. INC3 ACCIPN.MON' 

INCLUDE ' CAEGSDR ♦ INC J F LF.BLS . MON ' 

C 

DIMENSION X<4> ,X1<3> rX2<3> iY(4> rYl<3> * Y2<3> ?/AA) ? BUF(IOO) 
C 

DATA IREF/1/ 

C 

C**** CALCULATES THE AREA OF A GUADRALATERAL 
C 

IF ( NUNF’E • EG » 8 > GO TO 20 


C 

C 

C 


c 

r:; 

c 

20 


c 

30 


FOUR NODED ELEMENT 


XI 

(1) 

= 

X(1 

XI 

( 2 > 

= 

X ( 2 

XI 

(3) 

= 

X ( 4 

Y 1 

(1) 

= 

Y ( 1 

Y 1 

(2) 

= 

Y < 2 

Y 1 

(3) 

= 

Y ( 4 

A 1 

25 

EUTRIA 

X2 

(1) 

55 

X < 2 

X2 

< 2 > 

= 

X (3 

X2 

(3) 

55 

X ( 4 

Y2 

(1) 

= 

Y < 2 

Y2 

(2) 

= 

Y ( 3 

Y2 

(3) 

= 

Y ( 4 

A2 

55 

EUl 

RIA 


GO TO 30 


EIGHT NODED ELEMENT WITH STRAIGHT SIDES 


Xl(l) 

= 

xa 

X 1 ( 2 ) 

= 

X < 3 

XI ( 3 ) 


X < 5 

Yl(l) 

= 

Y ( 1 

Y 1 < 2 ) 

= 

Y ( 3 

Y1 (3) 

55 

Y ( 5 

A 1 = 

EUl 

RIA 

X 2 ( 1 ) 

= 

X<1 

X2 < 2 ) 

sr 

X ( 5 

X2 < 3 ) 

55 

X</ 

Y2 ( 1 ) 

= 

Y(1 

Y2 ( 2 ) 

55 

Y < 5 

Y 2 ( 3 ) 

25 

Y ( / 

A2 = 

EUT 

RIA 

EAG = 

A1+A2 


>U) 


,Y2> 



GO TO 900 


C 

807 

C 

877 

VOO 

C 


PRINT 877 * IEKR 

FORMATC 1X> 7 ACCELL* RETURNED WITH ERROR' r 14) 

CONTINUE 

RETURN 

END 
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SUBROUTINE COMP < PSI B, NT , NELM > 

CP**************** *********#******************************.***** 
CP# 

CP* comp: branches mj the appropriate ELEMENT TYPE Til CALCULATE 
CP* 

CP*************************************** ********************** 
CP* 

cp* description: 

CP* 

CP* ' COMP ' BRANCHES TO THE APPROPRIATE ELEMENT TYPE TO 

CP* CALCULATE THE COMPLIANCE AND THE COMPLIANCE SENSIT- 

CP* IVITY VECTOR ♦ 

CP* 

CP******************* **************************** ************** 

C 

INCLUDE ' C AEG SDR . INC 3 IMPL1C.SPC' 

INCLUDE ' C AEGSDR • 1NC!J ACC1PN. MON' 

INCLUDE ' CAEGSDR . INC3 CNTL . MON ' 

INCLUDE * CAEGSDR . INC! ELEDES , MON ' 

INCLUDE ' CAEGSDR « INC J SVECT K . MON ' 

COMMON/LCSDES/DLCS( 90 ) 

C 

EQUIVALENCE < NDAT (V7) f IDBS) , < NDAT < V8 > r I DBL ) 

C 

DIMENSION DAI N < 50 ) , PSIBB< bOO ) 9 PSIB ( 2 ) , CFBUh ( 6 ) ? PSIBTB ( 500 ) 
C 

DATA IREF/1/ 

C 

PSIB16 = 0*0 
PSIBC5 =0,0 
SDPSIT =0,0 
SDPSIB =0,0 
SDPSIH =0.0 
SDPSTB =0.0 
C 

C SET LOAD CASE NUMBER FOR COMPLIANCE 

r 

Ll = LCS + NC 

c 

0 

C SETUP POINTERS 
C 

CALL ACCELM ( 1 , IPNELM? IDBS r 1 ?0 r I ERR) 

IF < I ERR . NE . 0 ) CO TO 800 

CALL ACCFES (1*1 PNFfc S r I DBL r 1 r L 1 9 0 r 0 r IERR > 

IF(IERR.NE.O) GO TO 801 

CALL ACCCND (1*1 PNCNfl ? 1 DBL * 1 * Ll * 0 * 0 * I ERR ) 

IF < IERR . NE . 0 ) GO TO 802 

CALL ACCLCS< 1 * I PNLCS * I DBL * 1 » 0 * IERR ) 

IFUERR.NE.O) GO TO 805 

CALL ACCNOD ( 1 * J PNNOJJ 9 IDBS * 1 * Or IERR) 

IF ( IERR . NE . 0) GO TO 806 

CALL ACCELC< Ir 1FNELC, IDBS, 1,0,0, IERR) 

IF(IERR.NE.O) CO TO 807 

CALL ACCMAT <1,1 PNMAT r IDBS , 1 , 0 , 0 , I ERR ) 

IF ( IERR . NE ♦ 0 ) GO TO 808 

CALL ACCEPR ( 1 » I PNEPR r I DBS , 1 , 0 , 0 , J EKR ) 

IF ( IERR . NE . 0 > GO 10 809 

CALL ACCEENd r 1FNEEN r JDBL , 3 , L 1 , 0 , 0 , I ERR > 

IF< IERR.NE.O) GO TO 810 
C 
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C LOOP THROUGH THE ELEHENlS 
C 

DO 100 1 = 1 ? NELM 

IF( I ♦ UT • 1 > GO TO 50 

C 

C GET INTERNAL LOAD CASE NUMBER 
C 

CALL ACCLCS<2> IPNLCSf LI f2f DLCSr IERR) 

I F ( I ERR ♦ NE ♦ 0 ) GO TO 805 
ILCN = DLL'S ( 21 ) 

C 

C GET ELEMENT DESCRIPTORS 

C 

SO CALL ACCELh< 2 r JF’NELMf I f 2 f IEDf IERR > 

IF ( IERR # NL . 0 ) GO TO 800 
C 

C BRANCH TO THE APPROPRIATE ELEMENT TYPE 

C 

IF(ITYP.EQ.ll) CALL C0MP1 1 (N f f I r II. CN ) 

I F ( I T YP . EG « 5 ) CALL COMPOS < N T f I f ILCNf PSI BB ) 

IF ( ITYP * EG * 16 > CALL CP16(N T , I f ILCNf PSIBTB) 

C 

c 

PSIB16 = PSIB16+F*SIBTB ( I ) 

PS I DCS = PSIE<C5TPSIE<B< I ) 

SDPSIT = SDPSIT+riPS.IT(I) 

SDPSIF = SDPSIB+DPSIBCI) 

SDPSIH = SDPSIH+UPSIHU > 

SDPSTB = SDPSTBTDPSITB(I) 

100 CONTINUE 

<: 

IF (NT * GT * 1 ) GO TO 770 
WRITE < 1 0* 859 ) 

DO 710 1=1 f NELM 

710 WRITE < 10 f 857) I rDPSIT < I ) fDPSIB< 1 ) * DPSIH( I ) r DPS I TEC 1 > 
WR I TE ( 1 0 ? 861 ) SDPSIT » SDPSIBr SDPSIH f SDPST B 
C 

720 IF ( ITYP • EG « 1 1 ) CALL LPSI1 1 (PS1BT t ILCNfIFEF ) 

C 

PSIB(NT) = PSIBT + PSIBC5 + PSIB16 
C 

WRITE ( 10 • 858 ) PSIB(Nl) 

PRINT 858 r PSIB(N I) 

C 

C CLEAN-UP EVERYTHING 
C 

CALL ACC ELM (4 r) FNELN f O f 0 r 0 9 IERR > 

IF (IERR # NE* 0) GO TO 800 

CALL ACCFES ( A t I F’NF ES f 0 f 0 f 0 f 0 f 0 f 1 ERR ) 

IF ( IERR . NE * 0) GO TO 801 

CALL ACCCND ( A . 1PNCNU f 0 f Of Of 0 • 0 f IERR ) 

IF < IERR # NE . 0) GO TO 802 

CALL ACCLCS ( A t IFNLLS vOfOfOfl ERR ) 

I F ( I ERR « NE ♦ 0 ) GO TO 805 

CALL ACCN0D(4f IPNNOD 9 O 9 O 9 O 9 1ERK) 

IF < IERR . NE * 0) CO TO 806 

CALL ACCELC ( A , IPNELC f 0 f 0 f 0 f 0 f IERR ) 

IF ( IERR* NE * 0> GO TO 807 

CALL ACCMAT ( 4 r IPNMAT f 0 f 0 f 0 f 0 • I ERR ) 

IF(IERR.NE.O) GO TU 808 

CALL ACCEPR ( A t IPNEPR f 0 f 0 f 0 f 0 f I ERR > 
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IF ( I ERR ♦ Nt • 0 ) GO TO 80? 

CALL ACCEEN<4f IPNEEN.OfOf Of Of Of IERR) 

IF < IERR ♦ NE • 0 > GO TO 810 
C 

GO TO 820 
C 

C WRITE ERROR MESSAGES 10 THE SCREEN 
C 

800 PRINT 870 f IERR 

GO TO 820 

801 PRINT 871 r IERR 

GO TO 820 

802 PRINT 872 f IERR 

GO TO 820 

805 PRINT 875 f IERR 

GO TO 820 

806 PRINT 877 f IERR 

GO TO 820 

807 PRINT 878 f IERR 

GO TO 820 

808 PRINT 87? f IERR 

GO TO 820 

80? PRINT 876 f IERR 

GO TO 820 

310 PRINT 880 f IERR 

GO TO 820 

C 

c 

820 CONTINUE 

t; 

c 

857 FORMAT <I3f4Xf4<El6*8f4X) ) 

858 FORMAT C IX f ' PSIB* ' f E16 . 8 ) 

85? FORMAT UXf/UXf 'EN' ?6X? 'SENSITIVITY I ' f 7X f 'SENSITIVITY R'f 
*7Xf 'SENSITIVITY H ' f 6X f ' SENSI1 1 VI I Y I*') 

861 FORMAT UXf/>lXf'TOTAL=' f4(El6*8f4X> ) 

362 FORMATUXf 'ELEMENT 'f!4> 

970 FORMAT UXf'ACCELM RETURNED WliH ERROR '»14) 

871 FORMAT (lXf ' ACC FES RETURNED WITH ERROR 'fl4> 

872 FORMATUXf 'ACCCND RETURNED WITH ERROR XA> 

875 FORMAT < 1 X* ' ACCLCS RETURNED WITH ERROR ' t X4) 

876 FORMATUXf 'ACCEPR RETURNED WITH ERROR UX4) 

877 FORMATUXf 'ACCNOD RE fURNED WITH ERROR 'f*4> 

878 FORMATUXf 'ACCELC RETURNED WITH ERROR M4) 

87? FORMATUXf 'ACCMAT RETURNED WI lH ERROR 'fl4> 

880 FORMATUXf 'ACCEEN RETURNED WITH ERROR ' r 1 4 ) 

1004 FORMAT < 14 > 

C 

C 

c 


RETURN 

END 
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SUBROUTINE COMPOS <N (?I? ILCNjPSIBB) 

CP*************** ************ ********************* ************* 
CP* 

CP* compos: calculates compliance and sensitivity pop a flam 

CP* 

CP************************************************************* 

CP* 

cp* description: 

CP* 

CP* ' COMPOS ' CALCULATES THE COMPLIANCE ANH THE DESIGN 

CP* SENSITIVITY FOR A 1-D BEAM IN BENDING* WITH AN 

CP* APPLIED ELEMENT FORCE IN */IN. SELF WEIGHT IS 

CP* NEGLECTED. BEAM TORSION HAS BEEN ADDED* 

CP* 

CF ************************************************** *********** 
CP* 

CP* NT COUNTER FOR FINITE DIFFERENCE 

CP* I EXTERNAL ELEMENT NO. BEING PROCESSED 

CP* ILCN INTERNAL LOAD CASE NO. 

CP* NELM TOTAL NO. OF ELEMENTS 

CP* 

CP************************************************* ************ 


C 

INCLUDE ' CAEGSDR. J NC J IMPL1C.SPC' 

INCLUDE ' CAEGSDR. INC3 ACCIF’N.MON' 

INCLUDE ' CAEGSDR . INC J LNlL.MON' 

INCLUDE ' CAEGSDR. INC3 ELEDES.MON' 

INCLUDE ' CAEGSDR. INCT SVECTR . MON ' 

C 

EQUIVALENCE (ND AT (14) ? IP) 

C 

DIMENSION GF’LW ( 3 ) r PSI BB ( 500 ) > DDSHPF < 12 ) ? X < i ) r Y ( J > f 2 < 3 ) r 

* SHF'F (12)5 SBUF ( 200 ) * DAI N< 50 ) * CF BUF ( A ) t EQBUP ( 200 ) , 

* BUF ( 100 ) t CD ( 2 r 6 ) v U rU< i4 ) 9 C ( A r 2 ) v CDL ( 2 * 6 ) » T<3t 3) t 

* TB(6t 6) fC00R(3?3) 

C 

DATA GPLW/-.7745?66/D0f .ODOr .77459667P0/ 

DATA WTW/ . 55555556D0 ? .8888888VD0f .55555556110/ 

DATA KT/3/ 5 1 REF/1/ f MPT/1/ 

C 

PSIBB(I) = 0.0 
DF’SIBG = 0.0 
DF'SIHG = 0.0 
IF(I.GT.l) GO TO 50 
C 

C GET APPLIED FORCE IN */IN. 

C 

F*R I N T * f * 7 

PRINT *5 'ENTER APPLIED FORCE IN #/IN. UNITS' 

READ ( 5 5 1001 ) AF 
C 

C GET APPLIED DISTRIBUTED MOMENT 
C 

PRINT * 5 ' ' 

PRINT *f 'ENTER APPLIED DISTRIBUl ED MOMENI' 

READ ( 5 f 1 001 ) AM 

GET AREA MOMENT OF INERTIA ABOUT Y-AXIS 
50 CALL ACCEF'R < 2 , IPNEPR 5 1PTAB , 0 5 BUF 5 LF N* IERK* ) 
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IF< IERR. Nt ♦ 0) GO TO 80? 

YI = PUF ( 5 ) 

H = 2 # DOYBUF ( V ) 

P = 2 • DO*E<UF < 1 0 ) 

BW< N f ) = B 
BH(NT) = H 
C 

C GET WEIGHT DENSITY AND MODULUS OF ELASTICITY 
C 

GAMMA = 0.000 
E = PUF (5) 

0 * PUF ( 7 > 

G = E/<2.D0*<1 . DO+V ) ) 

C 

c 

C GET DISPLACEMENTS AT ELEMENT HMDS 


C 

DO 150 J= 1 > NUNF'E 

CALL ACCCND < 2 ? 1PNCND t IN f NN< J) » 1 y ILCN • CFBUF y LFN y J ERR ) 
IF < IERR. NE.O) GO TO 802 
DO 150 K* 1 9 NDOF 

CD(JfK) = CFBIJF < K ) 

150 CONTINUE 

C 

C********************************************************* * 


C EVALUATE DISFLS. AND CURVATURE AT THE GAUSS POINT USING 
C SHAPE FUNCTIONS - ONE FT. FOR CURV.r THREE PT. FOR DISPL 
C **************************************** ********* ********* 


C GET X, Yf AND Z OF ELEMENT NODES 

c; 

CALL ACCELC(2> 1PNEI CrKINl ? IREK ? PUF ? LEN t IERR ) 
IF ( IERR .NE.O) GO TO 807 
M a 1 

DO 200 J^ly9y3 
K = J+l 
L = J+2 
X< M ) = BUF(J) 

Y<M) * PUF(K) 

Z< M ) = PUF(L) 

M « h+1 

200 CONTINUE 

DO 210 J=ly3 

COOR ( 1 y J) “ X(J) 

COOR ( 2y J) = Y< J ) 

COOR < 3? J ) = Z(J> 

210 CONTINUE 

C 

C FORM THE ELEMENT LOCAL COORDINATE SYSTEM 
C 

IN3 = IN f NN ( 3 > 

CALL EUBTM < IN3 1 BE 1 A y COOR y T > IERR ) 

CALL ZEROSF* <TBf 34# IP) 

DO 220 J= 1 y 3 
DO 220 K~ 1 y 3 

TB(JyK) * 1 ( J f K > 

TP ( J+3 r K+3 > * T(JrK) 

220 CONTINUE 

CALL UMXAPTCTPf CD fCf 6y?y 6) 

DO 230 J= 1 f 2 
DO 230 K*1 >6 
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CDL(JtK) * C(KfJ) 

230 CONlINUE 

C 

C CALCULATE ELEMENT LENGTH 
C 

DX = X( 2 > -X < 1 ) 

DY = Y<2)-Y<1> 

DZ = Z< 2>-Z(l) 

EL * DSQRT < DX*DX+DY#DY+DZ#DZ ) 

If<r»Y*EO.O.ANrMfZ,EO*0) GO TO 236 
DO 235 J=1 1 NUNF’E 

CDL(Jf4) = -CPL(J>4) 

235 CONTINUE 
C 

C CHANGE LOCAL Y-ROTATION FROM POSII I ME TO NEGATIVE 
C IF PEAM LIES ALONG X GLOBAL AXIS 
C 

236 IF< DX * LT . 0* 001 ♦ AND « DX# GT « -0 • 001 ) GO TO 246 

DO 240 J=1 r NUNPE 

CDL ( J ? 5 ) = -CDL(JrS) 

240 CONI INUE 

C 

246 F = -AF - GAMMA*H*H 

C 

C CALCULATE THE TWISTING ANGLE 
C 

WXY = DABS< a:I.iL<2»4)-CDL(l?4> )/EL> 

C 

C EVALUATE SHAPE FUNCTIONS FOR IUSPL. - THREE POINT QUADRATURE 
C 


P2 

= B*B 

P3 

= B2*P 

P4 

= B3*B 

H2 

= H*H 

H3 

= H2*H 


WRITEUO* 1002) 1 
URITEUO* 1004) EfGjV 
WRITE<10t 1005) HfB 
C 

DO 300 K= 1 f 3 
PSI * GPLU(K) 

CALL EU3DSB < PSI > SHPF * DDSHPF > 2 * EL ) 

W = < SHPF ( 3 > #CUL (1*3 > +SHPF T 5 > #CDL (1*5) 

* +SHPF < 9 ) *CDL ( 2 f 3 > +SHPF <13 )*CI>L<2f 5> > 

C 

C EVALUATE SHAPE FUNCTIONS FOR CURV* - THREE POJNl CiUA C'RATU kE 
C 

WXX = (DDSHPF<3)#CDL<1* 3> +DDSHPK(b)*CDL< 1 » 5 > 

* + DDSHPF < 9 ) #CDL ( 2f 3 ) +DDSHPF ( 1 .1 >*CDL<2*5> ) 
WRITEdOf 860) KfUvUXXrUXr 

C 

IF(NT.GT.l) GO TO 250 
C 

C CALCULATE SENSITIVITY VECTORS 
C 

PJB - H3/3 « DO— 4 42P0¥B# ( H2+B4/ ( 4 « D0*H2 > ) 

PJH = B*H2-*42DO#E*2*<H-H4/d2.DOfH3) ) 

DPS IPG * DPSIBG+ ( -2*GAMMA#HYW- < E*H3/ 12 ) ♦UXXtfWxX- 

* P JB*G*WX Y*WX Y ) *U TW ( K ) * ( t L/2 . DO > 

DPSIHG » DPSIHG+(-2*GAhHA*B*W-<3*E*B*H;?/:l 2) *WXX#WXX- 

* PJH*G*WXY*WXr >*WfW<K)*<EL/2.P0) 



c 

C CALCULATE PSKB) - INTEGRAL OF FORCE*D I SPLACEMEN f 
C 

250 PSIBB ( I > = FSIBB(l) + <F*W+Ah*WXY >*WTW< K ) * < EL/2 . liO) 
300 CONTINUE 

WRITE< 10f 1003) DFSIBGf DPSIHG 
C 

IF ( NT #61*1) GO TO 820 
DPSIBO) = DPS I EG 
IipSIH(I) = DPSIHG 
C 

GO TO 820 
C 

C WRITE ERROR MESSAGES TD THE SCREEN 
C 


802 

PRINT 

872 9 

IERR 


GO TO 

820 


807 

PRINT 

878 f 

IERR 


GO TO 

820 


(308 

PRINT 

879 f 

IERR 


GO TO 

820 


809 

PRINT 

876 , 

IERR 


GO TO 

820 



c 

c: 

820 CONTINUE 

C 

C 

851 FORMAT (/f IX? 'BEAM WIDTH B=' 8. 5r 2X r ' BEAM L»LVTH« ' > F8 . 5> 2X 

*r ' E= ' * E9 « 3 > 2X > 7 1 Y r« 7 1 E9 . 3 t 2X f 7 GAMMA 35 ' rFA. 5» 7 APPLIED FORCE 
*=' fF8.5) 

855 FORMAT (IX? 'NODE*' , 12? 2Xt ' X*' f E12 * 5 1 2X > ' Y= ' * IK 12 . 5> 2X t ' 2- ' t 
*E12.5» 2Xf ' RX= ' f E 1 2 « 5 ? 2 X r 'KY* ' »£12«br2Xr 'RZ«' *fc 12.5) 

860 FORMAT < IX > ' GP* ' f I2f AX r ' W = ' t El 1 . 5 f 4X ? ' WXX= ' r El 1 ♦ 5 f 4X * 

* ' WXY= ' r El 1 • 5 ) 

070 FORMATUX* 'ACCELM RE) UR NED WITH ERROR 'tX4> 

872 FORMAT < 1 X ? ' ACCCND RE I UR NED WITH ERROR SH) 

876 FORMAT ( lXf 'ACCEPR RETURNED WITH ERROR S H) 

878 FORMAT < IX i ' ACLELC KE1 UKNED WITH ERROR '.J4> 

879 FORMAT (IX* 'ACCMAT Rt TURNED WITH ERROR 'rJ4> 

1001 FORMAT ( E12 ♦ 5 ) 

1002 FORMATdXf 'ELEMENT ='rt4> 

1003 FORMAT ( IX * ' L»PSI BG* ' ? E12 . 5 • 4X f ' DPSIHG* ' •FlI-'.S) 

1004 FORMAT ( IX f 'E* ' f El 2. 5r 4X r ' G= ' f E12 . 5 * 4X f ' V* ' rh.12 . 5> 

1005 FORMAT < IX * ' HEIGH 1 = ' i E 1 2 . 5 r 4X f 'WIDTH* ' f E12*5) 

C 

C 

C 

RETURN 

END 
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SUBROU r I NE COMP 11 < N f f I f ILCN > 

CPttt ******* ************************* ********************* ***** 
CH'fc 

CP# COMPUI CALCULATES COMPLIANCE AND SENSIT. FOR PL ANh STRESS 
CP* 

CP************************************************** *********** 
CP# 

cp# description: 

CP# 

CP# 'COMF1 1' CALCULATES THE COMPLIANCE AND INK DESIGN 

CP* SENSITIVITY OF THE FOUR AND EIGHT NODE PLANE STRESS 

CP* ELEMENT IN TRACTION WITHOUT SELF WEIGH f ♦ 

CP* 

CP***************************************************** ******** 
CP# 

CP* NT COUNTER FOR FINITE DIFFERENCE 

CP* I EXTERNAL ELFMENf NO. BEING PROCESSED 

CP* ILCN INTERNAL LOAD CASE NO. 

CP* 

CP************************************************************* 

C 

INCLUDE 'CAEGSDR. INCI IMPL1C.SPC' 

INCLUDE ' f AEGSDR * INCH ACCIPN.MON' 

INCLUDE ' C AEGSDR . INC3 CNTL.MON' 

INCLUDE 'C AEGSDR. 1NCJ ELFDES . MON ' 

INCLUDE ' C AEGSDR « INC3 SVtCTR.MON' 

C 

DIMENSION X<8> f Y<8) fZ<8> fSHPFCH) fGPL<2f4) ,DATN(50> f 

* BUF < 1 00 ) f PSI B < ?) * SBUh ( bO > fCFBUF <6> fEOBUT ( t>0 > r 

* DSHPGX(B) fDSHPGY(B) fBF <4f 40> r SE< S oO > f DSHPL O' f H > f 

* SIGMA ( 6 f 4 ) f EPSLN ( 6 f 4 ) 

C 

DATA GPL/2*- . 57735027 f . 5773502/ f - . 5/7 35027 f 

* 2* . 57735027 f - . 57735027 f . 57735027/ 

DATA KT/3/ f IREF/ l/ 

C 

SE < I > * 0. 

IF ( NT ♦ GT ♦ 1 ) GO TO 350 
C 

C GET ELEMENT STRESSES AND STRAINS 
C 

CALL ACCFES ( 2 • IPNTESf KIN I f I REF • ILCN f SBUF f I ENf I ERR ) 

IF ( I ERR . NE . 0 ) GOTO HOI 
LOC * LEN - 1 
M = 1 

DO 50 K=1fNSVAL 

SIGMA(IfK) = SRUF(M) 

SIGMA < 2 f K ) = SBUF (M+l ) 

SIGMA ( 3 f K > = SBUF (M+3) 

M = M+4 

50 CONTINUE 

DO 60 K“1fNSVAL 

EPSLN(IfK) = SBUF < M) 

EF'SLN ( 2 f K ) = SBUF (M+l) 

EPSLN ( 3 f K ) = SBUF (M+3 ) 

M = M + 4 

60 CONTINUE 

C 

C GET X AND Y FOR JACOBIAN EVALUATION 
C 

CALL ACCELC(2f IPNELCfKINTf IREFfBUFfI ENBf JFRR) 
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IF < I ERR . NE * 0) GO TO 807 
M = 1 

DO 200 L = 1 1 LENBf 3 
J = L+l 
K = L+2 
X < H ) = BUF < L ) 

Y(H) = BUF(J) 

Z(h> - BUF <K) 
h = h + 1 

200 CONTINUE 

C 

C 

C CALCULATE FORCES AT THE GAUSS POINTS 
C 

DO 250 L= 1 f NI'OF 
DO 250 K= 1 r NSVAL 
BF < K f L ) = 0.0 
250 CONTINUE 

C 

C LOOP OVER THE GAUSS POINTS 
C 

DO 300 K~1 f NSVAL 
PSI = GPL(lfK) 

ETA - GPL(2vK) 

C 

C EVALUATE SHAPE FUNCTIONS AT THE GAUSS POINTS 
C 

IF ( ISTYP « EQ # 2 > LAI L EU2DI. Q ( PSI • L i‘A . K I * SHPF ? DSHPL f 

* DSHPGX 9 DSHPGY f DET J ? X f Y f I ERR) 

IF < ISTYP. EG « 4) CALL EU2DPQ (PSIf E LA . K f ? SHPF , DSHKL f 

* DSHPGX 7 DSHL'GY f »E TJ f X f Y r I ERR > 
IFCIERR.NE.O) GOTO 809 

300 CONTINUE 

WRI TE ( 1 0 f 855 ) 

DO 320 K=*1f NSVAL 

320 WRI TE ( 1 0? 854 ) Kf < SIGMA < J f K > f J= 1 f NSIG > 

WRITE ( 1 0 f 860 ) 

DO 330 K=1f NSVAL 

330 WRITE ( 1 0 f 854 ) Kf ( L PSLN < J f K > f J= 1 f NS I G ) 

DO 340 J=1fNSIG 
DO 340 K» If NSVAL 

SE < I ) = SE( I > + SIGHA(J'k)*EPSLN(»l9K>*D£TJ 
340 CONTINUE 

C CALCULATE SENSITIVITY VECTOR 
C 

DPS IT ( I ) = - SEC I) 

350 GO TO 820 
C 

C WRITE ERROR MESSAGES TO THE SCREEN 
C 


801 

PR] 

[NT 

871 f 

I ERR 


GO 

TO 

820 


007 

PR] 

[NT 

878 f 

I ERR 


GO 

TO 

820 


809 

PR] 

[NT 

Q76f 

I ERR 


GO 

TO 

820 



C 

C 

820 

C 


CONTINUE 



n n o 


c 

854 FORMAT < IX , 12 , 2X , 3 (E16 . 8 » i!X ) ) 

855 FORMAT ( IX, 'GP' ,5Xr 'SIGMAX ( OP) ' > 8X» 'SIGMAYXGP) ' »8X, 
# ' SIGMAXY ( GP ) ' ) 

860 FORMAT < IX, ' GP' » 5X, ' EPSLNX'GP) ' ,8X , 'KPSLNY (OP) ' , 8X , 

*'EPSLNXY(GP> ' ) 

871 FORMAT ( IX, ' ACCFES REJURNED WITH ERROR ',14) 

876 FORMAT < IX, 'EU2HPQ RETURNED WITH ERROR ',14) 

878 FORMAT < IX, ' ACCELC RETURNED WITH ERROR ' ,X4> 


RETURN 

END 
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SUBROUTINE CP16 < N T> I r ILCN* PSIBTB ) 

CP ************************************************************* 
CFt 

CP* CP16 5 BRANCHES TO THE APPROPRIATE ELK MEN I SUBTYPE 
CP* 

CP************************************************************* 

CP* 

CP* DESCRIPTION t 
CP* 

CP* 'CPI A' BRANCHES TO THE APPROPRIATE ELK MEN f SUBTYPE 

CP* TO CALCULATE THE COMPLIANCE AND COMPLIANCE HESTON 

CP* SENSITIVITY OF THE PLATE BENDING ELK MK N T 16. 

CP* NOTE: THIS DOES NOT TAKE INTO ACCOUNT ANY MEMBRANE 

CP* STIFFNESS. 

CP* 

CP************************************************************* 

CP* 

CP* NT COUNTER FOR FINITE DIFFERENCE 

CP* I EXTERNAL ELEMENT NO. BEING PROCESSED 

CP* ILCN INTERNAL LOAD CASE NO. 

CP* PSIBTB COMPLIANCE 01" A BENDING PLATE t USED FOR 

CP* CALCULATING THE FINITE DIFFERENCE 

CP* 

CP* **************************************************** ******** 
INCLUDE ' CAEGSDR. INC3 IMPLIC.SPC' 

INCLUDE 'CAEGSDR. INC J ELIDES. MON' 

INCLUDE ' TAEGSDR. I NCI SVECTR.MON' 

C 

DIMENSION PSIBTB(SOO) 

C 

C** BRANCH TO THE APPROPRI ATK ELEMENT SUBTYPE 
C 

IF(ISTYP.EQ.l) CALL CPI 601 ( N T r I r ILCN ? PSIBTB ) 

IF ( I ST YP . EQ . 2 > CALL CPI 602 ( N T v I r ILCN > PS 1 B TB ) 

C 

RETURN 

END 



SUBROUTINE CPI 601 (NT r I f ILCN r PSIBTB ) 

CP ************************************************************ 
CP* 

CP* CP1601J CALCULATES COMPLIANCE AND DESIGN SENSITIVITY 
CP* 

CP************************************************** ********** 

CP* 

cp* description: 

CP* 

CP* ' CP 1601 ' CALCULATES COMPLIANCE AND DESIGN SENSITIVITY 

CP* FOR A TRIANGULAR BENDING FLhMENT • 

CP* 

CP************************************************************ 


CP* 

CP* 

CP* 

CP* 

CP* 


NT COUNTER FOR FINITE DIFFERENCE 

I EXTERNAL. El t MEN T NO. BEING PROCESSED 

ILCN INTERNAL LOAD CASE NO. 

PSIBTB COMPLIANCE f USED FOR CALCULATING FINITE DIFFERENCE 


CP* 

CP********************* *********************************** ***** 

c 

INCLUDE / C AEGSDR. INC3 IMPL1C.SPC' 

INCLUDE ' C AEG SDR * INC 3 ACOIPN.MON' 

INCLUDE ' C AEGSDR .INCH CNlL.MON' 

INCLUDE 7 C AEGSDR . INC J ELEDES . MON * 

INCLUDE ' C AEGSDR .INC.! SVECTR.MON' 

C 

DIMENSION XCO f Y( 3 ) f PSIB< 2 ) f SBUF <100 ) * CPE ( 500 > fEF<J.6> f 

* SIGMA (6> fEPSLN( 6) r BUT ( 1 00 ) f CFBUM 6 ) ? EQBUF < 100 ) * 

* PSIBTB <500 rCD(3»6> fZC5) 

C 

DATA KT/J/ f IREF/J /r MPT/J / 

C 

CPE < I ) s 0. 

PSIBTB ( I ) = 0. 

C 

C GET PROPERTIES 

C 

CALL ACCEPR(2? IPNEPRf IPTABfOf BUF fLENf IERR ) 

IF ( I ERR ♦ NE . 0) GO TO 80? 

PB ( NT ) = BUFC>5) 

IF ( NT « GT. 1 ) GO TO 65 
C 

C GET ELEMENT STRESSES AND STRAINS 

C 

CALL ACCFES < 2 r IPNFES r KINT f ) REF r ] LCN r SBUF f Lfe N f I ERR ) 

IF < I ERR * Nb . 0 ) GOTO 801 
M = 1 

DO 50 J=1fND0F 

SIGMA ( J > = SBUF(J) 

50 CONTINUE 

M = 7 

DO 60 J=1fND0F 

EPSLN(J) = SBUF ( M > 

M = M-f 1 

60 CONTINUE 

C 

C GET DISPLACEMENTS FOR PSI CALCULATION 
C 

65 DO 70 J= 1 f NUNP t. 

CALL ACCCND ( 2 f IP'NCNDf INTNN ( J > f 1 f ILCN f CK HUT fIENf IERR) 



IF ( I ERR * NE ♦ 0) GO TO 302 
DO 70 K*1 9 NDOF 

CD(JfK) = CFBUF(K) 

70 CONTINUE 

C 

C GET EGUI VALEN f FORCES AT THE ELEMENT NODES FOR PS1 CALC • 
C 

CALL ACCEEN< 2 1 1PNEEN * KIN I ? I REF r 1LCN > EOBUF ? LEN • IERR) 
IF ( IERR * NE ♦ 0 > GO TO 808 
M = 1 

DO 80 J*1»NUNPE 
DO 80 K=1 f NDOF 

EF(JrK) = EGBUK(M) 

M * M + l 

30 CONTINUE 

IF<NT,GT,1> GO TO 350 
C 

C GET THE JACOBIAN 
C 

CALL ACCELC ( 2 1 IF'NELC r KIN f * IREF * BUF t LFNB f IERR ) 

IF ( I ERR . NE . 0 > GO TO 807 
h = 1 

DO 200 J-1,LENFU3 
X (M ) = BIJF(J) 

Y ( M ) = BUFCJ+l) 

Z ( M ) = BUF < J+2) 
h = M+l 

200 CONTINUE 

DETJ = EUTRIA(XrV) 

C 

C 

C CALCULATE SENSITIVITY VECTOR 
C 

DO 340 J=1 f 3 

CFE(I) = CF'E(I> + SIGMA( J)*EPSI.N< J)*DETJ 
340 CONTINUE 

DPSITB(I) = - CPE < I ) 

C 

C CALCULATE F'SI(B) - INTEGRAL OF FORCE*DI SPI ALFMF N T IN l 
C: 

350 DO 400 J=1 t NIJNF’E 

FSIBTBU) = PS IBUtKI) + EF( Jt 3>*CD( J,3) 

400 CONTINUE 
GO TO 820 
C 

C WRITE ERROR MESSAGES TO THE SCREEN 
C 


801 

PRINT 

871 f 

IERR 


GO TO 

820 


002 

PRINT 

872 > 

IERR 


GO TO 

820 


807 

PRINT 

877 1 

IERR 


GO TO 

820 


008 

PRINT 

879 t 

IERR 


GO TO 

820 


80? 

PRINT 

880 1 

IERR 


GO TO 

820 



C 

C 

820 

C 


CONTINUE 



cs n 


853 

871 

872 
877 
87? 
880 
t: 

C 

c 


FORMAT < 1 X f 'NODE' f6Xf 'DISP X' fOXf ' DISP 
* j 8X f ' ROT X ' f 9X ? ' POT Y'fVXf'RUT 

FORMAT <1XfJ4f6(2XfE12*4>) 

FORMAT < IX f ' ACCFES RETURNED WITH ERROR 
FORMAT < 1 X f ' AC CCND RETURNED WITH ERROR 
FORMAT (IX f 'ACCELC RE I UKNED WITH ERROR 
FORMAT ( IX f ' ACCEEN RETUKNUD WITH ERROR 
FORMAT < IX f ' ACCEPK RETURNED WITH ERROR 


Y' fbXf 'DISP 

/' ) 

' r :m ) 

' f X 4 ) 

' r M> 

' f34) 

' f M ) 


RETURN 

END 
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SUBROUTINE CP160?<Nf * ] t ILCN? PSIBTB ) 

CF ■*****************************************+****************** 

CF'* 

CP* CR1602J CALCULATES COMPLIANCE ANN DESIGN SENSITIVITY 
CP* 

CP******************************************* ***************** 

CP* 

cp* description: 

CP* 

CP* 'CP 1602 ' CALCULATES COMPLIANCE AND DESIGN SENSITIVITY 

CP* FOR A FOUR NODE BENDING ELEMENT. 

CP* 

CP************************************************************ 

CP* 

CP* Nf COUNTER FOR FINITE DIFFERENCE 

CP* I EXTERNAL NO. BEING PROCESSED 

CP* ILCN INTERNAL LOAD CASE NO. 

CP* PSIBTB COMPLIANCE 

CP* 

CP************** **************************************** ******* 

C 

INCLUDE ' CAEGSEiR . INCH IMPLIC.SPC' 

INCLUDE ' C AEGSDR ♦ INC3 ACCIPN.MON' 

INCLUDE * CAEGSDR. INC J C N 1 L . MON ' 

INCLUDE 'i AEGSDR.INCD ELEDES. MON ' 

INCLUDE ' CAEGSDR. INC J SVECTR.MON' 

C 

DIMENSION X<4> ? Y(4> rPSIB (2) , SBUF (1 00) ,CPF <500 > > EF ( A 1 6 ) jI < A ) , 

* SIGMA < 6 f 4 ) f EF’SLN ( 6 ? 4 > *BUF”(JOO)*CFBUf (6) * 

* EQBUE ( 3 00 > r PS I BTB ( 500 ) > CD < * > 6 > 

C 

DATA KT/3/f IREF/l/r MPT/1/ 

C 

CPE < I ) = 0. 

PSIBTB(Z) = 0. 

C 

C GET PROPET I ES 

C 

CALL ACCEPR < 2 > 1 FNLPR > IPTAB f 0 f BUF > LEN r I ERR ) 

IF( I ERR. N£ • 0) DO TO HO? 

F'B ( NT > = BUF (25) 

IF ( NT . GT . 1 > GO TO 65 
C 

C GET ELEMENT STRESSES AND STRAINS 

C 

CALL ACCFES( ?» IPWKES* KIN! r I REF 9 I LCN f SBUF" r I EN r I ERR) 

IF < I ERR # NE . 0 ) CO 10 801 
M = 1 

DO 50 K~1 f NSVAL 
J = M+l 
L = M+2 

SIGMA ( 1 f K) * SBUF (M) 

SIGMA ( 2 f K > = SBUK<J> 

SIGMA ( 3 f K ) * SBUF” < L ) 

M « M + 6 

50 CONTINUE 

DO 60 K*l» NSVAL 
J = H+l 
L = M+2 

EPSLN(lrK) = SBUF ( M ) 



EPSLN< 2 f K ) = SDUF(J) 

EPSLN < 3 * K ) = SBUF CL) 
h = M+6 

60 CONTINUE 

C 

C GET DISPLACEMENTS FOR PSI CALCULATION 
C 

65 DO 70 JolrNUNPfc" 

CALL ACCCND(2f IPNCNDy INTNN( J) f 1 * 1LCN» CFBUK * LF N* J ERR) 
IF ( I ERR « NE ♦ 0 > t=0 TO 802 

DO 70 K=1 9 NOOF 

CD(JfK) = CFHUK(K) 

70 CONTINUE 

C 

C GET EQUIVALENT FORCES AT THE ELEMENT NODES FOR PSI CALC ♦ 

C 

CALL ACCEEN ( 2 * IPNEFN* KIN I 1 1 REF * ILCN? EQBUh »I.ENr 1ERR) 

IF ( I ERR « NE ♦ 0) CO TO 808 
M « 1 

DO 80 J-1»NUNPE 
DO 80 K“1 f NDOF 

EFCJfK) = EGBUF <M> 
h = M + l 

80 CONTINUE 

I F ( NT * GT * 1 ) GO TO 350 

C 

C GET THE JACOBIAN 

c: 

CALL ACC EL C ( 2 t I F'NELC • K IN f * 1 RFK r BUF * l.ENB > I ERR ) 

IF( I ERR « NE * 0) CO TO 807 
M = 1 

DO 200 J-l ? LENH » 3 
X < M ) = BUF < J> 

Y ( M ) = BUF ( J+ 1 ) 

2 < M ) = BUF<J+2> 

M = M+l 

200 CONTINUE 

DETJ = AREAQ < X f Y > 

DETJ = DETJ/ 4 * DO 
C 
C 

C CALCULATE SENSITIVITY VECTOR 
C 

DO 340 J=1 t 3 

DO 340 K=1 f NSVAL 

CPE < I ) * CPE(I) + SIGMA < J f K ) *EPSI.N < J» K ) YDET J 
340 CONTINUE 

DPSITBCI) = - CPE ( I ) 

C 

C CALCULATE PSI(B) ~ INTEGRAL OF FORCE*D I SPL ACFMEN T IN Z 
C 

350 DO 400 J= 1 * NUNPE 

PSIDTB(I) = PS 1 BIB < I ) + t>‘< J'3>*CD< J*3> 

400 CONTINUE 
GO TO 820 
C 

C WRITE ERROR MESSAGES TO THE SCREEN 
C 

801 PRINT 871* I ERR 

GO TO 820 

002 PRINT 872 * I ERR 
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*307 

808 

*309 

C 

c 

820 

C 

C 

852 

85i 

*371 

872 

877 

879 

880 
C 

C 

C 


GO TO 

820 


PRINT 

877 * 

I ERR 

GO TO 

820 


PRINT 

879* 

I ERR 

GO TO 

820 


PRINT 

880* 

I ERR 

GO TO 

820 


CONTINUE 


FORMAT ( IXf 

' NODE 

t 

*8X 

* 'KOT 


* *8X*'K0T X ' f 9X * ' ROT Y'fVXf'ROT /') 
FORMAT <1X*I4*6<?X*E12*4) ) 

F ORMA r<lXf ' ACCFES KEIURNED WITH ERROR ' *14) 
FORMAT (IX*' ACCCND RETURNED WITH FRROR '*14) 
FORMATOX* 'ACCELC RETUKNtD WITH ERROR '*14) 
FORMAT ( IX r 'ACCEEN RE I UR NED WITH ERROR '*J4> 
FORMAT < IX? ' ACCEPR RETURNED WITH ERROR '*14) 


RETURN 

END 



SUBROU f INE CPS 11 1 ( PSIBT H 1. 1 » J REF > 

C 

INCLUDE 'CAEGSDR. JMCJ INPLIC.SPC' 

INCLUDE ' C AEGSDR .INC! ACCIPN.HON' 

INCLUDE ' CAEGSDR# INC 3 CNTL . WON ' 

INCLUDE 7 TAEGSDR* INC3 KLEDES ♦ MON ' 

C 

DIMENSION DATN ( 50 ) f CFBUF < 6 ) 

C 

C CALCULATE F SKFO FOR PLATE - FORCFSDISPL . ACEMh N I 
C 

PRINT *?' ' 

PRINT * f ' ENTER NODE WHERE LOAD IS APPLIED ' 

READ< 5r 1 004 > NEXT 
PRINT *f' ' 

PRINT # f ' ENTER LOAD DIRECTION (1 2 X f • . ♦ f 6 J K*7 ) ' 

READ < 5 f 1004 ) LMIR 
C 

C GET INTERNAL NODE NUMBER 
C 

CALL ACCNOD < 2 r IPNNOD f NEX1 » 2 f DATN f I ERR ) 

IF< IERR.NE.O) GO TO 806 
NINT = D A TN ( 4 ) 

C 

C GET DISPLACEMENT AT NODE 
C 

CALL ACCCND < 2 ? .1 PNCND f NINT f IREFf IL 1 f OFDUH . I. EN . 3 FRR ) 
IF (IERR.NE.O) CO TO 802 
DISP = CFBUF (LDIR) 

C 

PSIBT = 40000*DISP 

C 

C WRITE ERROR MESSAGES 
C 


802 

PRINT 872 f I ERR 
GO TO 820 




1306 

PRINT 877 f I ERR 
GO TO 820 




C 





020 

C 

CONTINUE 




c 

072 

FORMAT (IXf 'ACCCNJ i 

RETURNED 

WITH 

ERROR ' f j 4 ) 

877 

F0RMAT(1Xf 'ACCNOD 

RET URNKl) 

UI1H 

ERROR ' f X 4 ) 

1004 

C 

FORMAT (14) 




C 

RETURN 

END 






SUBROUT I N E HI SR ( PSIB > NT r Nt.LM ) 

CP************************************************************ 

Cp# 

cp# pisp: branches to the appropriate element type for the 
cp* 

CP*#******#*# ************************************ ************* 
CP* 

cp* description: 

CF‘* 

CP* 7 DISP ' BRANCHES TO THE* APPROPRIATE ELEMENT TYRE FOR 

CP* CALCULATION OF THE DISPLACEMENT CONSTRAINT AND THE 

CP* DISPLACEMENT SENSITIVITY VECTOR, 

CP* 

CP* ************************* ******************** ************** 
CP* 

CP* PS IB DISPLACEMENT AT X 

CP* NT COUNTER FOR FINITF DIFFERENCE ~ 

CP* NELM TOTAL NO, OF ELEMENTS 

CP* 

CP************************************************************ 

C 

INCLUDE ' CAEGSDR * INC 3 IMPLIC.SPC' 

INCLUDE ' CAEGSDR, INCH ACCIPN.MON' 

INCLUDE ' f AEGSDR, INC 3 CN fL , MON ' 

INCLUDE ' CAEGSDR, INC 3 ELFPKS, MON' 

INCLUDE ' CAEGSDR, INC3 SVECTR.MON' 

C 

EQUIVALENCE < NDAT <97 ) ? I DBS ) ? (NHAT ( 98) f J DDL ) 

C 

DIMENSION DA TN< 50 > f PSI B < 2 > f CFBUF (6) 

C 

DATA IREF/1/ 

C 

SDPSTB =0.0 
SUPS IT =0,0 
SDPSIB = 0,0 
SDPS I H =0,0 
IF (NT , GT , 1 ) DO TP 10 
C 
C 

C SET ADJOINT LOAD CASE NUMBER 
C 

L2 = LCS + NC 
C 

C SET ORIGINAL LOAD CASE NUMBER 
C 

10 LI = 1 
C 

c; 

C SETUP POINTERS 
C 

CALL ACCELM ( 1 t IPNELMf IDBS f 1 f Of IERR) 

IF < IERR . NE , 0 > GO TO 800 

CALL ACCFES ( 1 r I PNKFSf I DBL f 1 f LI f 0 f 0 ? IERR ) 

IF(IERR.NE.O) GO TO 802 

CALL ACCCND < 1 1 I PNCNH f I DBL f 1 f L 1 f 0 r 0 f IERR ) 

IF( I ERR , ME , 0 ) GO TO 802 

CALL ACCNOD ( 1 f I PNNOB f I DBS r 1 f 0 f I ERR ) 

IF(IERR.NE.O) GO TO 806 

CALL ACCELC ( 1 f I PNELC f J DBS f 1 » 0 r 0 f I ERR ) 

IF(IERR.NF.O) GO TO 807 



84 


CALL ACCMAT ( 1 r IPNMAT > I DPS 7 1 > 0 7 0 > ) EKR > 

IF< IERR* NE . 0) GO TO 008 

CALL ACCEPR ( 1 f IPNEPRf I DBS? 1 vOf OvIERR) 

IF ( IERR * NE « 0) GO TO 909 

c: 

C LOOP THROUGH THE BUFFERS TO OFT STRESSES r STRAINS ? hOMEN fS ? ETC « 
l 

DO 100 1=1 ? NELM 

c 

C GET ELEMENT DESCRIPTORS 
C 

CALL ACCELM<2* lPNEI.hr If JREFr IFDr IERR) 

I F ( I ERR » NE ♦ 0 ) GO TO 800 
C 

I F ( NT . GT « 1 ) GO TO 720 
C 

C BRANCH TO THE APPROPRIATE EL If MEN ‘I TYPfc 
C 

IF ( I TYP • EG* 1 1 ) CALL DISP11 <N I > I r LI f LV' )‘L J ) 

IF ( I T YP * EQ « 5 ) CALL DISP05(Nl 7 I 7 LI tL2f )L1) 

IF < I TYP ♦ EG* 16) CALL DP16 <N T r I r LI * L2* IL1 ) 

C 

C 

SDPSIT = SDPSIT+HPSIT(I) 

SDPSIB * SBPSIB+riPSIBd) 

SUPS I H = SDPSIH+DPSIH ( I ) 

SDPSIB = SDPSTBTDPSI IB < I ) 

100 CONTINUE 
C 

WRITE ( 10 j 859 ) 

DO 710 1 = 1 f NELh 

710 WRITE (10* 857) J t DPSI T ( I ) * DPSIB ( I ) f DI'SIH ( J ) rlipSITB(I) 

WRITE < 1 0 7 86 1 ) SDPSIT 7 SDPSI B 7 SDPSIH • SDPSI B 
C 
t 

C CALCULATE PSI(B) - ORIGINAL DISPL AT NODE UHERF AD J I* I.UAD APPLIED 
C 

720 PRINT tf ' ' 

PRINT *, ' ENTER NODE NUMBER WHERE ADJl LOAD IS APPLIED' 

READ (5? 1004 ) NEXT 
LDIR = 3 
C 

C GET INTERNAL NODE NUMBER 

C 

CALL ACCNOD <2 7 IPNNOD t NEXT r DATNr IFRR) 

IF ( IERR * NE ♦ 0) GO TO 806 
HINT = DATN< 4 ) 

C 

C GET DISPLACEMENT AT NODE 
C 

CALL ACCCND < 2 7 IPNCND » NIN I r IREF r IL1 r CFBUF f LFN r IFRR ) 

IF < IERR ♦ NE . 0) GO TO 802 
DIS = CFBUF (LD I R ) 

c 

PS IB ( NT ) = -DIS 

c 

WRITE (iOt 858) PSIB(Nf) 

PRINT 850 1 PSIB(NT) 

C 

IF (NT • EG. 1 ) GO TO 730 
WRITE ( 107*) 



WRITE (10? 856) Nt XT 


C 

C 

C CLEAN-UP EVERYTHING 
C 

730 CALL ACCELM <4,3 PNEI M, 0,0, 0,1 ERR ) 

IF < I ERR * NE • 0 > GO TO 800 

CALL ACCFES <4,1 PNFES ,0,0, 0,0, 0,1 ERR ) 

IF < I ERR « Nb ♦ 0 > GO TO 801 

CALL ACCCND <4,1 PNC NO , 0 , 0 , 0 , 0 , 0 , 1 ERR ) 

IF < 1ERR « NE • 0) GO TO 802 

IF (NT . GT * 1) GO TO 750 

CALL ACCLCS <4,1 PNLCS ,0,0,0, 1 ERR ) 

IF < I ERR • NE « 0 > GO TO 805 
750 CALL ACCNOD <4,1 PNNOO ,0,0,0, 1ERR ) 

IF < I ERR . NE * 0) GO TO 806 

CALL ACCELC < 4 , 1 PNEI C , 0 , 0 , 0 , 0 , 1 ERR ) 

IF < I ERR * NE * 0) GO TO 807 

CALL ACCMAT <4,1 PNMAT , 0 , 0 , 0 , 0 , 1 ERR ) 

IF < I ERR ♦ NE ♦ 0) GO TO 808 

CALL ACCEPR <4,3 PNEPR , 0 , 0 , 0 , 0 , 1 ERR ) 

IF < I ERR * NE * 0 > GO TO 809 
C 

GO TO 820 
C 

C WRITE ERROR MESSAGES TO THE SCREEN 
C 


800 

PRINT 

870, 

I ERR 


GO TO 

820 


801 

PRINT 

871, 

IERR 


GO TO 

820 


802 

PRINT 

872, 

I ERR 


GO TO 

S20 


805 

PR IN I 

875, 

IERR 


GO TO 

820 


806 

PRINT 

877, 

IERR 


GO TO 

820 


807 

PRINT 

878, 

IERR 


GO TO 

820 


808 

PRINT 

879, 

IERR 


GO TO 

820 


80? 

PRINT 

876, 

IERR 


GO TO 

820 



C 

c 

820 CONTINUE 


C 

c 

H56 FORMAT < IX , ' #**ADJ01N f LOAD IS APPLIED AT NODE', 14) 

85/ FORMAT < I3,4X,4<E16*8,4X) ) 

858 FORMAT < IX , ' PSIB« ' ,E16 * 8 ) 

85? FORMAT < IX,/, IX, 'EN' ,6X, 'SENSITIVITY T' , 7X , ' SENSITIVITY H ' 

*7X, 'SENSITIVITY H',6X, 'SENSTIVITY TH' > 

861 FORMAT < IX,/, IX, ' TOTAL- ' , 4 < E 16 . 8, 4X ) ) 

862 FORMAT < IX, 'ELEMEN f ',14) 

870 FORMAT < IX, 'ACCELM RETURNED WITH ERROR ',34) 

871 FORMAT < IX, 'ACCFES R L TURNED WITH ERROR ',)4> 

872 FORMATUX, 'ACCCNU RETURNED WITH ERROR ',14) 

875 FORMAT < IX, ' ACCLCS RETURNED WITH ERROR ',34) 

876 FORMAT < IX , ' ACCEPR RETURNED WITH ERROR ',34) 

877 FURMAT< IX, 'ACCNOD RETURNED WITH ERROR ',14) 
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!378 FORMAT < IX » ' AOCELC RETURNED WITH ERROR ',14) 

879 FORMAT < IX , 'Al;CM AT RETURNED WITH ERROR ',14) 

1004 FORMAT (14) 

2000 FORMAT ( A ) 

C 

C 

C 


RETURN 

END 
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SUBROUTINE DISP05<N f f I r LI . I 2* 11. 1) 
CP************************************************************* 
CP* 

CP* DISF’05! CALCULATES THE DISF’L. DESIGN SENSITIVITY Or A BEAM 
CP* 

CP*********** ************************ ************************** 
CP* 

cp* description: 

CP* 

CP* 7 DISP05 ' CALCULATES THE DISPLACEMENT CONSTRAINT DESIGN 

CP* SENSITIVITY OF A 1-P BEAM IN BENDING i WITH AN APPLIFD 

CP* ELEMENT FORCE IN */JN. SFLF WEIGHT )S NEGLECTED . 

CP* TORSION IS ACTIVE* 

CP* 

CF’*** **************************** *************** *************** 


CP* 




CP* 

NT 

COUNTER 

FOR FINITE DIFFERENCE 

CP* 

I 

EXTERNAL 

ELEMENT NO. PEING PROCESSED 

CP* 

LI 

ORIGINAL 

EXTERNAL LOAD CASE NO. 

CP* 

L2 

EXTERNAL 

ADJOINT L.OAD CASE NO. 

CP* 

CP* 

IL1 

INTERNAL 

LOAD CASE NO. OF ORIGINAL L.OAD 


CP************************************************************ 

C 

INCLUDE ' CAEGSDR . INC3 IMPLIC.SF'C' 

INCLUDE ' CAEGSDR. INC3 ACCIPN.MON' 

INCLUDE ' CAEGSDR. INCI CNTL . MON 7 
INCLUDE 7 CAEGSDR. INCH ELEDES* MON 7 
INCLUDE 7 CAEGSDR. INCH SVECTR.MON' 

COMMON/LCSDES/LiLCS<?0> 

C 

EQUIVALENCE (NDAT ( ] 4 > r IF > v ( NDAT ( VS > r I DDL ) 

DIMENSION X<3> * Y<3> fZ. < 3) * DDSHpF (1 2 ) > DATN ( 50 ) . BUK (I 00 ) * 

* SHPF( 12) fCFBUF (A) .GF*LW<3> . Al D( 2 > 6 ) > CD ( 2 >6 ) . W fW ( 3 > r 

* C < 6 . 2 ) . P < 6 . 2 > . T < 3 . 3 ) ? TB < 6 ? 6 ) .LDL(2<6> r ALDL (2.6) . 

* COOR (3.3) 

C 

DATA GPLW/-. 7/459667 r .0 0000000- .77459667/ 

DATA WTU/ .555555556. .80888889. .55555556/ 

DATA KT/3/f IREF/l/f MPT/3/ 

C 

DPSIBG =0.0 
DPSIHG =0.0 
IF(I.GT.l) GO TO 50 

c; 

C REQUEST APPLIED FORCE IN LOAD/LENGTH 
C 

PRINT *.' ' 

PRINT * . ' ENTER APPLIED LOAD IN C FORCE /l. ENGTHI 7 
READ< 5. 1001 ) AF 

C 

C GET AREA MOMENT 01" INERTIA ABOUT Y-AXIS 
C 

50 CALL ACCEPR < 2 . IPNEF’R .IF‘TAB.0. BUK . LEN. IERR) 

IF (IERR.NE.O) GO TO 809 
YI = BUF ( 5 ) 

H = 2* BUK (9) 

B = 2* BUF < 10 > 

PH (NT) = H 
BW ( N I ) = B 
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c get weight density and modulus of elasticity 
c 

CALL ACCMAT < 2 * IPNMAT r NMAT ? MPT f BUF f LEN f I ERR ) 

IF < I ERR • NE ♦ 0 ) GO TO 808 
GAMMA = 0 « ODO 
E = BUF < 5 > 

V = BUF (7) 

G = E/(2«D0*(1 . DO+V) ) 

C 

IF ( I • 6T • 1 ) GO TO 60 
C 

C GET INTERNAL LOAD CASE NUMBER FOR ORIGINAt LOAD 
f; 

60 CALL ACCLCS(IfIPNLCSfIDBLfIfOfIERR) 

IF ( IERR.NE.O) GO TO 805 

CALL ACCLCS ( 2 f IPNLCS f L 1 f 2 f DLCS f I ERR ) 

IF ( I ERR* ♦ NE » 0 ) GO TO 805 
IL1 = DLCS ( 21 > 

C 

C GET DISPLACE MEN I S AT ELF MEN f ENDS 
C 

CALL ACCCNO ( 1 f IPNCNDf IDBL f 1 f LI f 0 f Of I ERR > 

IF(IERR.NE.O) CO TO 802 
DO 70 J=1fNUNPE 

CALL ACCCND<2f IPNCNDf INTNN( J) f 1 f IL1 rCFDUh fL.ENfIFRR) 
IF ( I ERR ♦ NE ♦ 0 ) GO TO 802 
DO 70 K=1fNU0F 

CD(J-K) ■ CFBUF(K) 

70 CONTINUE 

C 

t: GET INTERNAL LOAD CASE NUMBER FOR ADJOINT LOAD 

C 

CALL ACCLCS < 1 f IPNLCSf IDBL f 1 fOf IERFO 

IF ( I ERR . ME ♦ 0 ) GO TO 805 

CALL ACCLCS<2f IPNLCS f L2 f 2 f DLLS f IERR) 

IF ( I ERR * NE . 0 ) GO TO 805 
ILCN = DLCS (21) 

C 

C GET DISFLACEHCN fS AT ELEMENT ENDS FOR ADJOINT LUAD 
C 

CALL ACCCND (1 f I PNC Nil , ] DBL f!fL2fOfOfIERR) 

IF( IERR « NE ♦ 0 ) GO TO 802 
DO 100 J=1fNUNPE 

CALL ACCCND ( 2 f IPNCNDf IN TNN( J) f 1 f /LCNfCFDUFfLENf IERR) 
IF(IERR.NE.O) GO TO 802 
DO 100 K«1fND0F 

ALD(JfK) « CFBUF(K) 

100 CONTINUE 

c 

C*************** *********************************** ******** 

C EVALUATE PISPLS. AND CURVATURE AT THE GAUSS POINT USING 
C SHAPE FUNCTIONS - ONE PT. FOR CUKV* f 1'HkEE PT. FOR IMSPL 
C* ********** ******** ***** ********************************** 

C GET Xf Yf AND Z OF ELEMENT NODES 
C 

CALL ACCELC(2f IPNELCf KINT f IREFfBUKfLENf IERR) 

IF( IERR « NE « 0> GO TO 807 
M = 1 

DO 200 J=1 *9f3 
K = JO 
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X ( H ) = BUF<J> 

Y<M) = BUMK) 

Z(H) = BUF(L) 

M = Mil 

200 CONTINUE 

DO 210 J=1 f 3 

COOR(lfJ) * X(J) 

COOR ( 2 f J ) = Y(J) 

COOR ( 3 f J ) « Z(J) 

210 CONTINUE 

C 

C FORM THE ELEMENT LOCAL COORDINATE SYSTEM FOR D1SPLACEMEN I S 
L 

IN3 = IN I NN (3 ) 

CALL EUBTM ( I N3 r BETA • COUR f I f 1 ERR > 

CALL ZEROSPC TBf36*IP> 

DO 220 J=1 » 3 
DO 220 K=1 f 3 

TB< Jf K) a 1 ( J f K ) 

TB < J43 * Ki3> « T(J*K> 

220 CONTINUE 

CALL !JMXABT( TBtCDf C?6?2j6) 

DO 230 J= 1 f 2 
DO 230 K»1 f6 

CDL(JfK) = C(Kf J> 

230 CONTINUE 

CALL UMXABT ( I B f ALD • U f 6 ? 2 f 6 > 

DO 232 J=1 f 2 
DO 232 K«1 ?6 

ALDL(JfK) = D(KfJ) 

232 CONTINUE 

C 

C CALCULATE ELFhEN \ LENGTH 
C 

DX = X(2>~X(1) 

DY = Y(2)-Y(l> 

DZ = Z(2)-Z(l> 

EL = DSQRT<DX*DX+DY*DY4DZ*DZ> 

C 

C CHANGE LOCAL Y-ROTA I ION FROM POSITIVE TO NEGATIVE .1 F 
C BEAM LIES ALONG X GLOBAL AXIS 
C 

IFCDX.LT.O. 001 .AND. DX.GT. -0.001 > GO TO 246 
DO 240 J= 1 f NUNh'E 

CDL(Jf5> = -CDL<Jf5> 

ALDLCJfS) = -ALDL < J f 5 ) 

240 CONTINUE 

C 

246 F = -AF - GAMMA*B*H 

C 

C CALCULATE THE TWISTING ANPLFS 
C 

UXY a DABSC <CDL<2 f4)-CDL<1f4> >/EL) 

AWXY a DABS ( ( ALDL (2*4) -ALDL (lr4))/bL) 

C 

C EVALUATE SHAPE FUNCTIONS FOR PISFL* - THREE PUINf QUADRATURE 
C 

B2 a b*B 
B3 = P2*B 
B4 a B3*B 
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H 2 * H*H 
H3 = H2*H 
C 

DO 300 K=lr3 
F'SI = UP LU(K) 

CALL EU3DSB(PSIfSHPFf DDSHPF* 2f EL) 

W « <SHPF<3)*CDL< 1 * 3 > +SHPF< 5) *CDL ( 1 r5)+SHPF<9>* 

* CDL < 2 * 3 ) +SHPF (11) *CDL ( 2 r 5 > ) 

AW = (SHPF ( 3 ) #ALDL (If 3 >+SHPF ( 5) ftALDL (If t») +SHPF ( 9 ) * 

* ALDL(2f3)+SHPF< 1 1 >*ALDL< 2* 5 > ) 

C 

C EVALUATE SHAPE FUNCTIONS FOK OUUV. - I HKEE POINT QUADRATURE 
C 

WXX a (DDSHPF(3)*CDL<ly3)+)irSHPF(5)#CHL(l»5)+ 

* DDSHPF < 9 ) *CDL ( 2r 3 ) +DDSHPF ( J 1 ) *CDL < 2 * 5 > > 

AWXX = < DDSHPF < 3 )*ALDL(1 r 31 +DDSHPK<5) *AI.DL< 1 * 5> + 

* DDSHPF ( 9 > *ALDl <2 f 3 ) +UDSHPF (11) *ALDL < 2 f 5 > > 

C 

C CALCULATE SENSITIVITY VECTORS 
C 

PJB = H3/3 « DO- *42D0*B*(H2+B4/(4 ♦ t»0*H2 ) ) 

PJH = B*H2-.42DO*B2*(H-M/<12.DO*H3> ) 

DPS I BG « IjPSIBG+<-GAMMA*H*AW-<E*H3/12>*AWXX*UXX- 

* PJB*G*WXY*AWXY>*WrW<K>*(fr'L/2.B0) 

DPS IMG =* DPSIHG+<-GAMMA*B*AW-<3*E*B*H2/J2>*AWXX*WXX- 

* FJH*G*UXY*AUXY>*Wl W(K)* (EI./2.H0) 

300 CONTINUE 

DPSIBU) a DPSIBG 
DF'SI H < I ) « DPSIHG 
C 

GO TO 820 
C 

C WRITE ERROR MESSAGES TO THE SCREEN 
C 


302 

PRINT 

872 f 

IERR 


GO TO 

820 


805 

PRINT 

875 ? 

I ERR 


GO TO 

820 


807 

PRINT 

878 r 

IERR 


GO TO 

820 


808 

PRINT 

879 t 

IERR 


GO TO 

820 


(309 

PRINT 

876 f 

IERR 


GO TO 

820 



c 

c 


820 CONTINUE 

C 

C 

851 FORMAT ( / ? IX * 7 BEAM WIDTH B* ' f F8 « 5f 2X f ' BEAM DLPTH* ' f F8 « 5* 2X 

*r 'E*' fE9. 3?2X» 'IYY*' *E9.3 f 2X r 'GAMMA* ' fF6.5»2Xf ' APFLIED 
♦FORCE*' fF8.5> 

855 FORMAT ( IX f 'NODE*' » I 2r 2X » 'X*' 'E12'5r 2 X» ' Y*' fE12.5v2Xf 'Z*' f 
»E12.5f2Xf 'RX*' fE12.5f2Xf 'RY* ' f fcl2 .5 f2Xf 'RZ* ' fE12,5> 

860 FORMAT (IXf 'GP*' f I2f 4Xf ' W*' fUI . 5f 4Xf ' WXX* ' fEI 1 . 5f 4Xf ' AW* ' f 

♦ El 1 .5* 4X? ' AWXX= ' fEII .5) 

864 FORMAT ( IXf ' NODE- ' f)2?2Xf' AX*' f£'12#5f2Xf 'AY* ' rE12. 5>2 Xf 
*'AZ=' fE12.5>2Xf 'ARX-' fE12.Sf2Xf 'ARY-' *E12,5f2Xf ' A R Z « ' * 

* E12.5) 

1372 FORMAT (IXf ' ACCCND RETURNED WITH ERROR 'f!4> 

875 FORMATUXf 'ALCLCS RETURNED WITH ERROR '»I4) 


C-3- 



(376 FORMAT < IXf ' ACCERR RETURNED WITH 

879 FORMAT < 1 X j ' ACCELC RETURNED WITH 

87? FORMAT < 1 X * ' ACCHAT RETURNED WITH 
1001 FORMAT (El 2 # 5> 

C 

c 

RETURN 

END 


error ',ia> 

ERROR '?XA> 
ERROR ' tXA) 
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SUBROUTINE DISP1 1 <Nf> I »L1 »l.2» IL1) 
cm*********************************************************** 
CF* 

CP# IHSPlli CLACULATES THE DISFL* GENSIT* FOR PLANE STRESS 
CP# 

CP* **************************** ######*##♦#*##**##### *********** 
CP# 

c:p* description: 


cp 

CP 

(;P 

CP 

CP 


' DISF'l 1 ' CAL.CULATLS THE DXSRLACEMFNl CONSTRAINT DESIGN 
SENSITIVITY FOR THE FOUR AND EIGHT NODE PLANE STRESS 
ELEMENT WITH TRACTION ? SELF WEIGHT NOT INCLUDED* 


CP************************************************************* 

COUNTER FUR FINITE DIFFERENCE 
EXTERNAL ELEMENT NO* WING PROCESSED 
ORIGINAL EXTERNAL LOAD CASE NO* 

EXTERNAL ADJOINT l UAD CASE NO* 

ORIGINAL INTERNAL LOAD CASE NO* RETURNED TO DISP.EOR 
CP#**#***** ******************************************** ******** 

c 


CP* 


CP* 

NT 

CP* 

I 

CP* 

LI 

CP* 

L2 

CP* 

CP* 

I LI. 


INCLUDE 

INCLUDE 

INCLUDE 

INCLUDE 

INCLUDE 


* CAEGSDR . INC 3 
7 CAEGSDR# INC J 
' CAEGSDR . I NCI 
' CAEGSDR. INC J 
' CAEGSDR * INC 3 


IMPLIC* SP’C * 
ACCJPN# MON ' 
CNTL ♦ MON ' 
ELEDES * MON ' 
SVLXTR . MON * 


COMMON/LCSUES/DLCS < 90 ) 

EQUIVALENCE (HflAT(98> r 1DDL > 

DIMENSION SHPF<8) ?GPL(2t 4 ) t DU* ( 100) fX< 8> r Y<8> r Z<8> » 

* DATN(SO) rPSJB<2) fSBUK(SO) r USHPGX <8 ) r DSHFGY (8) > 

* DSHPLC'f 8) ? SIGMA (6 f A) t EPSLN< 6 f A > r DF ( 4 ,48 > > SE ( 500) 

DATA GPL/2*- . 5773502 7 » * 5/735027 > - * 5//35027 , 

* 2* * 57735027 1 - . 57735027 , * 5 7/35027/ 

DATA KT/3/ , I REF/ 1/ ? ME 7/1/ 


C 

c 

c 

10 


JJ » LI 

GET INTERNAL LOAD CASE NUNXt.N 

CAL L ACCLCS < 1 , J PNI .CS f I DDL ,1,0,1 ERR ) 

IF < I ERR * NE . 0 ) GO TO 905 

CALL ACCLCS < 2 1 IFNLCSf JJ ,2 , 1*1 CS , 1 ERR ) 

IF< I ERR . NE * 0) GO TO 805 
ILCN « DLCS< 21 ) 

SETUP POINTER FOR STRESS-STRAIN DUFFER 

CALL ACCFES ( 1 r I PNFES f ) DDL , 1 , ILCN ,0,0,1 ERR > 
IF ( I ERR ♦ NE * 0 ) GO TO 801 
IF ( JJ • EG « LI ) 1L1 * ILCN 

SE ( I ) = 0*0 


C 

C GET ELEMENT STRESSES AND STRAINS 
C 

CALL ACCFES ( 2 , IPNFES , KIN f , IREF , ILCN , SDUF , t EN , I ERR) 
IF( IERR.NE.O) GOTO 801 
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LOC * L EN - 1 
M * 1 

IF <JJ. EQ.L2) 60 TO 55 
I»0 50 K~1 fNSVAL 

SIGMA ( 1 f K ) * SBUF (M > 

SIGMA < 2 f K) * SBUK <M+1 ) 

SIGMA < 3 f K> = SBUF<M+3> 

M = M FA 

50 CONTINUE 

GO TO 65 
55 M = 17 

DO 60 K = 1 9 NSVAL 

EPSLN < 1 f K ) = SBUF(M) 

EPSLN ( 2 f K ) = SBUF(Mil) 

EPSLN ( 3 f K > - SBUF ( Mi 3 ) 

M = Mi 4 

60 CONTINUE 

65 IF <jj.eq.l:>> GO TO 100 

JJ = L2 
GO TO 10 
C 

C GET X AND Y FOR JACOBIAN EVALUATION 
C 

1 00 CALL ACCELC ( 2 f IPNELC f K IN f f I REF f HUH f l.ENB f I ERR ) 

IF ( IERR. NL'.O ) 60 TO 807 

M = 1 

DO 200 L=1fLENBf3 
X ( M ) = BIJF <L) 

Y(M) ^ BUF (L+l ) 

Z ( M ) = BUF (L+2 ) 

M = Mil 

200 CONTINUE 

C 

C CALCULATE FORCES AT THE GAUSS POINTS 
C 

DO 250 L^JfNDOF 
DO 250 K^IfNSMAL 
BF(KfL) = 0.0 
250 CONTINUE 

C 

C LOOP OVER THE GAUSS POINTS 
C 

DO 300 K*! f NSVAL 
PSI = GPL (It K) 

ETA * GPL(2fK) 

C 

C EVALUATE SHAPE FUNCTIONS AT THE GAUSS POINTS 
C 

IF ( ISTYP ♦ EG. 2) CALL EU2DLQ< PSI f ETA f K f f SHPH f DSHPL f 

♦ DSHPGX t DS HP GY f DET J f X f Y f 3 ERR > 

IF< ISTYP.EG. 4) CALL EU2DPQ< PSI f ETA r K T f SHPF f DSHPL f 

* DSHPGX f DSHPG Y f DFTJ f X f Y r IERR > 

IF ( IERR * NL . 0) GOTO HOV 

300 CONTINUE 

DO 340 J^IfNSIG 
DO 340 K~ If NSVAL 

SE( I ) ~ SE<I) + SIGMA( JfK)*FPSI.TM JfK)*D£TJ 
340 CONTINUE 

C 

C CALCULATE SENSITIVITY VECTOR 
C 
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DPS IT ( I > = - SE ( I > 

700 CONTINUE 
C 

GO TO 820 

C WRITE ERROR MESSAGES TO THE SCREEN 
C 

801 PRINT 871 f IERR 

GO TO 820 

805 PRINT 875 r IERR 

GO TO 820 

007 PRINT 878 f IERR 

GO TO 820 

009 PRINT 876* IERR 
GO TO 820 
C 
C 

820 CONTINUE 

C 

C 

854 FORMAT < IX* I2*2X>3<E16.8*2X) > 

855 FORMAT < IX* 'GP' »5Xr 'SIGMAX<GF*) ' * HX * ' SIGMAY ( GP ) ' *8X* 
♦ ' SIGMAXY ( GP > ' ) 

860 FORMAT < IX * 'GP' *5X* 'EPSLNX(GP) ' *8X* 'EPSLNY <(>P) ' * 8X* 

* ' EPSLNXY ( GP ) ' > 

871 FORMAT ( IX * ' ACCFFS RETURNED WITH ERROR SI 4) 

875 FORMAT (IX? 7 ACCLCS RETURNED UI(H ERROR '*I4> 

876 FORMAT < IX* ' EU2DP0 RETURNED WITH ERROR 'rX4) 

878 FORMAT ( IX * 'ACCELC RETURNED WITH ERROR '*14> 

l; 

C 

RETURN 

END 


SUBROUTINE DF*16 ( N f ? ) ,L 1. , L2> IL1 > 

CF •*#*******************************************************'**** 
CF'* 

CF'* CP16J BRANCHES TO THE APPROPRIATE ELEMENT SUBTYPE 
CP* 

CF : '*************** ********************************************** 
CP* 

cp* description: 

CF'* 

CP* ' CF‘16 ' BRANCHES TO THE APPROPRIATE ELEMENT SUBTYPE 

CP* TO CALCULATE THE DISPLACE HEN I DESIGN SENSITIVITY OF 

CF‘* THE PLATE BENDING ELEMENT 16* 

CP* note: this does not take into ACCOUNT amt membrane 

CP* STIFFNESS. 

CP* 

CF’* ************************************************************ 


CP* 
t p* 

NT 

CP* 

I 

CP* 

LI 

CP* 

L2 

CP* 

IL1 

CP* 



COUNTER FOR FINITE DIFFERENCE 
EXTERNAL ELEMENT NO* BEING PROCESSED 
ORIGINAL EXTERNAL LOAD CASE NO. 

EXTERNAL ADJOINT LOAD CASE NO. 

ORIGINAL INTERNAL LOAD CASE NO. - RETURNED VALUE 


CP******** ********************************************* ******** 


INCLUDE 'CAEGSDR.INC3 IhPLIC.SPC' 
INCLUDE ' C AEG SDR * INCH KL EDE3 . MON ' 
INCLUDE 7 C AEG SDR .INCH SVtCTR.MON' 


C 

C * * BRANCH TO THE APPROPRIATE ELEMENT SUBTYPE 
C 

IF < ISTYP ♦ EQ . 1 ) CALL DPI 601 (N f f I r LI ? L2f IL1 ) 
IF < ISTYP.EQ.2) CALL DP1602 ( NT» I »L1 r L2f 1L.1 > 

c 

RETURN 

END 


SUBROUTINE DPI AO 1 ( N f » ] rl.i » l.2» 1L1 ) 

CPt******************* ******************* ********************** 
CF* 

OF* DF1601J CALCULATES THE DISPLACEMENT CONSTRAINT DESIGN 
CF'* 

CP***** **************** ************** ************************** 

VF* 

cp# description: 

CP* 

CP* " DPI 60 1 / CALCULATES THE DJ SPL.ACEMENT CONSTRAINT DESIGN 

CP* SENSITIVITY VECTOR FOR A TRIANGULAR PLATE BEN 10 NO 

CP* ELEMENT « 

CP* 

CP ************************************************************* 
CP* 

CP* NT COUNTER FOR FINITE DIFFERENCE 

CP* I EXTERNAL ELEMENT NO. BEING PROCESSED 

CP* LI ORIGINAL. EXTERNAL LOAD CASE NO. 

CP* L2 EXTERNAL ADJOIN f* LOAD CASE NO . 

CP* IL1 ORIGINAL INTERNAL LOAD CASE NO. - RETURNED VALUE 

CP* 

CP************************************************************* 
1NCLVVE ' LAEGSDR . INC3 IMPLIC.SPC' 

INCLUDE ' C AEGSDR • INC J ACCIPN.hON' 

INCLUDE ' CAEGSDR ♦ INC! CNTL . MON ' 

INCLUDE ' CAEGSDR . INC3 ELFDES ♦ MON ' 

INCLUDE ' CAEGSDR ♦ INC J SVECT R . MON ' 

COHMON/L.CSDES/DLCS ( 90 ) 

C 

EQUIVALENCE ( NDAT ( 98) f I DDL > 

C 

DIMENSION SIGMA ( 6 ) f EPSLN ( 6 > t SBUh (100 > f »UF < x 00 ) f CPE < 500 ) f 
* X(3> rY(3> f2(3> 

C 

DATA KT/3/f IREF/1/ 

C 

JJ = LI 
CPE ( I ) =0.0 

C 

C*** GET INTERNAL LOAD CASE NUMBER 
C 

10 CALL ACCLCSO fIPNLCSfIDBLfIfOfIERR) 

IF(IERR.NE.O) GO TO B05 

CALL ACCLCS < 2 f IPNLCS fJJf2fDLCSf IERR) 

IF ( I ERR ♦ NE . 0 ) GO TO 805 
ILCN = HLCS(21) 

c 

C*** GET PROPERTIES 
C 

CALL ACCEPR < 2 f IPNEPR f IP TAP f 0 f BUF r LFN f J ERR > 

IF< IERR.NE.O) GO TO B09 
PB(NT) = BUF < 25) 

C 

C*** SETUP POINTER FOR STRESS-STRAIN BUFFER 
C 

CALL ACCFES ( 1 f IPNFES f IDBL f i ? X LCNr 0 f 0 f ) ERR ) 

IF ( IERR . NF . 0 ) 00 TO 801 

IF(JJ.EQ.Ll) I L 1 = ILCN 

C 

C*** GET ELEMFNl STRESSES AND STRAINS 
C 



50 

55 


60 

65 


C 

c*** 

c 

100 


200 

c 

c 

c*** 

c 


340 

C 

C 

c*** 

c 

801 
805 
GO 7 
G09 
C 

820 

t; 

854 

855 

860 

871 

875 


CALL ACCFES(2»XFNhFSrKJNT» IRFFf ILCNtSBUFf LENrlERR) 
ZF(IERR.NE.O) GOTO HOI 
IF ( JJ ♦ EG . L2 ) 60 TO 55 

DO 50 K*l»NJtKiF 

SIGMA ( K ) * SBUF(K) 

CONTINUE 
GO TO 65 
M « 7 

DO 60 KMfNDOF 

EFSLN(K) = SfeUK(h) 

M » M+l 
CONTINUE 

IF ( JJ*£CK L2> GO TO 100 
JJ = L2 
GO TO 10 

GET THE JACOBIAN 

CALL ACCELC <2* IPNELCf KINT ? IREK* HUF f LEND f I ERR ) 

IF( I ERR * NE » 0 ) GO TO 807 
h = 1 

BO 200 JalrLENR'3 
K * J+l 
LL = J+2 
X ( M ) = BUF(J) 

Y< M ) = BUF(K) 

2(H) = BUF(LL) 

M = M+l 
CONI INUE 

DETJ a EUTRIA<XrY> 

CALCULATE SENSITIVITY VECTOR 
DO 340 J=1 f 3 

CPE ( I ) * 1:PE<1> + SIGMA( J)*EFSLN< J)*D£TJ 
CONTINUE 

nPSITB(I) a -CPE ( I ) 

GO TO 820 

WRITE ERROR MESSAGES III THE SCREEN 

PRINT 871 f IERR 
GO TO 820 
PRINT 875 ? IERR 
GO TO 820 
PRINT 878 9 IERR 
GO TO 82 0 
PRINT 879 f IERR 
GO TO 820 

CONTINUE 


FORMAT (1XfI2f2Xf3(E16*8f2X> ) 

FORMAT ( IX f 'GP' f 5X» ' SIGMAX ( GP ) ' fHXf 'SIGMAY(GP) ' f8Xf 

* ' SI GMAXY ( GF‘ ) ' ) 

FORMAT < 1 Xf 'OF' f5Xf 'EFSLNX(GF) ' f8Xf ' EPSLNY ( GP ) ' fHXf 

* ' EPSLNXY < GP > ' ) 

FORMAT ( IX f 'ACCFES KF l UliNKu WITH ERROR 'fI4> 

FORMAT ( IXf / ACCLCS RETURNED WITH ERROR 'fI4> 



n n 


07 8 FORMAT ( IX, 'ACCELC RETURNED WITH ERROR ' t\A) 

87? FORMAT ( IX t 'ACCEF'R RETURNED WITH ERROR ',IA> 


RETURN 

END 
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SUBROUTINE DPI 602 (N T r I > U. r I. 2r 111 ) 
CP**************************************** ********************* 
CP* 

C F : ‘* DP 1602* CALCULATES THE DISPLACEMENT CONST RAIN I DESIGN 
CP* 

CP************************************************************* 

CP* 

cp* description: 

CP* 

CP* ' DP 1 602 7 CALCULATES THE DISPLACEMENT CONSTRAINT DESIGN 

CP* SENSITIVITY VECTOR FOR A POUR NODE PLATE BENDING 

CP* ELEMENT. 

CP* 

CP************************************************************* 


CP* 




CP* 

NT 

COUNTER 

FOR FINITE DIFFERENCE 

CP* 

I 

EXTERNAL 

ELEMENT NO. BEING PRUCFSSfc.il 

CP* 

LI 

ORIGINAL 

EXTERNAL LOAD CASE NO. 

CP* 

L2 

EXTERNAL 

ADJOINT LOAD CASE MO. 

CP* 

CP* 

I L. 1 

ORIGINAL 

INlERNAL LOAD CASE NO. - RETURNED 


CF ******************************************************* ******* 
INCLUDE ' CAEGSDR . INC J IMPLIC.SPC' 

INCLUDE ' CAEGSDR .INC! ACCIPN.MON' 

INCLUDE ' CAEGSDR. INC3 CNfL.MON' 

INCLUDE ' CAEGSDR. INC3 ELEDFS. MON ' 

INCLUDE ' CAEGSDR. INC J SVECTR. MON ' 

COMMON/LCSDES/DLCS< ?0 ) 

C 

EQUIVALENCE (NDAT(?8> r IDBL) 

C 

DIMENSION SIGMA ( 6r 4 > t EPSLN (6*4) f SBUT ( 100.) • BUF < 100 ) , CPE < 500 > f 
* X<4) , Y(4) fZ<4> 

C 

DATA KT/^/f IREF/1/ 

C 

JJ = LI 
CPE ( I ) =0.0 
C 

C*** GET PROPERTIES 

C 

CALL ACCEFR<2r JFNEPRr IPTAB r 0 > BUF r LFN f IFRR ) 

IP ( I ERR • NE . 0 ) GO TO SO? 

PB(NT) = BUF ( 25) 

C 

C*** GET INTERNAL LOAD CASE NUMBER 
C 

10 CALL ACCLCSC IfIFNI.CSf IDBL >1*0* I ERR) 

IF ( I ERR . Nt . 0 ) GO TO 805 

CALL ACCLCS ( 2 ? IPNLCSr JJ • 2 f UL.CS t I ERR > 

IF< IERR.NE.O) GO TU 805 
I LCN = DLCS ( 21 ) 

C 

C*** SETUP POINTER FOR STRESS-STRAIN BUFFER 
C 

CALL ACCFES< 1 r IFNFES f ItiBL* 1 f I L CN »0?0* TERR) 

IF (IERR.NE.O) GO TO 801 
IF(JJ.EQ.Ll) IL1 * I LCN 
C 

C*** GET ELEMENT STRESSES AND STRAINS 
C 
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CALL ACCFES (2*1 PNFES r K IN I , ] REF» ll.L’Nf SBUh » LEN • I ERR ) 
IF < I ERR . NE • 0 ) GOTO 601 
M = 1 

IF (JJ.EG.L?) GO TO 55 
DO 50 K*1 r NSVAI . 

J = M + l 
L = M+2 

SIGMA ( 1 f K ) = $><UK<h> 

S IGMA ( 2 > K > = SBUFCJ) 

SIGMACJfK) = SBUFCL) 

h = h + 6 

50 CONTINUE 

GO TO 65 
55 M = 25 

DO 60 K*1 f NSVAL 
J = M+l 
L = M+2 

EF’SLN < 1 > K > = SBUF(h) 

EPSLN < 2 > K ) = SBUF < J > 

EPSLN(3»K> = SBUFCL) 

M = M + 6 

60 CONTINUE 

65 IF <JJ.EQ.L2) GO TO 100 

JJ = L2 
GO TO 10 
C 

C *** GET THE JACOBIAN 
C 

100 CALL ACCELC < 2* XPNELC t KINT t IREK? BUF ? LENH » j FKR ) 

IF< I ERR . Nt . 0) GO TO 807 
M = 1 

DO 200 J = 1 r LENB r 3 
K = J+l 
LL = J+2 
X ( M ) « BUFCJ) 

Y(h) = BUF (K) 

Z < M ) = BUF ( LL) 

M = M + l 

200 CONTINUE 

C 

DETJ » AREAQ < X * Y > 

DETJ » DETJ/ 4 ♦ DO 

c: 

C*** CALCULATE SENSITIVITY VECTOR 
C 

DO 340 J=1 f 3 

DO 340 K=1 i NSVAL 

CPE ( I ) = CPE ( I ) + SIGMA ( J ? K ) *EFSLN < J t K ) *DET J 
340 CONTINUE 

DPSITB(I) = -CPE ( I ) 

C 

GO TO 820 
C 

C*** WRITE ERROR MESSAGES TO THE SCREEN 
C 

801 PRINT 871 ? I ERR 

GO TO 820 

805 PRINI 875 , I ERR 

GO TO 820 
PRINT 878 f I ERR 
GO TO 820 


807 
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809 

e 

820 

c 

854 

855 

{360 

871 

875 

878 

879 
C 

C 


PRINT 879 ? I E RR 
60 TO 820 

CONTINUE. 


FORMAT OXrI2r2Xr3(E16*8r2X) ) 
FORMAT ( IX f ' GP * r 5X •' SI GMAX ( GP ) 
* ' SIGMAXY ( GP ) ' ) 

FORMAT (IX r 'GP' ,5X* ' fc'PSLNX ( 6P ) 
% ' EPSLNXY ( GP ) f ) 

FORMAKlXr 'ACCFES RETURNED 
FORMAT < IXr 'ACCLCS RETURNED 
FORMAT ( IX r ' ACCELL KEIUKNED 
FORMAT < 1 X r ' ACCEFR RE fUNNED 


f8Xf 'SIGMAY (GP) ' rSX* 


? 8X f 'EPSLNY (UP) 7 ?BXt 


WITH 

h HR UR 

' * JA> 

WITH 

ERROR 

' ,XA) 

WITH 

ERROR 

' • .»* A ) 

WITH 

ERROR 

' ,1 A) 


RETURN 

END 


SUBROUTINE EU^I'SB ( RSI > SURF > UDSHKF t KT » EL ) 

CF****** ***************** ***************** *******'* ************* 
CF'# 

CP* EU3DSB5 BEAM SHAPE FUNCTIONS 
CP# 

CP# ********************************* *************************** 


CP* 



CP* 

PSI 

GAUSS POINT 1. OCA’I ION ♦ THE CENTER OF THE BF AM 

CP* 


IS ZERO. 

CP* 

SHPF 

STANDARD BEAM SHAPE FUNCTIONS - RETURNED VAI UL 

CP* 

DDSHPF 

SECOND DERIVATIVE (JF THE SHAPE FUNCTIONS 

CP* 


- RETURNED VAl.Uli 

CP* 

KT 

FLAG FOR RETURNING DDSHPF 

CP* 


=1 ONLY RETURN SHPF 

CP* 


=2 RETURN BOTH SHPF AND DDSHPF 

CP* 

CP* 

EL 

ELEMb'Nl LENGTH 


CP******************************************* ****************** 
C 

INCLUDE ' C AEGSDR . INC 3 IMPLIC.S PC' 

C 

DIMENSION SHPF ( 12 ) t DDSHPF (12) 

C 

EL2 = EL*EL 
EL3 = EL2*EL. 

X « PSI + (EL/2) 

X2 = X*X 
X3 « X2*X 
t; 

C CALCULATE THE SHAPE FUNCTIONS 
C 

SHPF(l) = t-X/EL 

SHFF<2> ■ 1-(3*X;VEL2> M2*X3/EI 3> 

SHPF (3 > = SHKF( 2) 

SHPF < 4 ) * SHPF U ) 

SHPF ( 5 ) = X-(2*X2/EL) KX3/EL2) 

SHPF (6) = SHpF (5) 

SHPF ( 7 ) = X/t"L. 

SHPF ( 8 > = < 3 * X2 / EL V > -- ( 2 * X 3/ EL 3 ) 

SHPF ( 9 > = SHPF(S) 

SHPF (10) « SHPF (7) 

SHPF (11) = (-X2/EL ) i (X3/EL2 ) 

SHPF (12) = SHPF (11) 

C 

IF (KT . EG . 1 ) 00 TO 900 

C 

C CALCULATE THE SECOND DERIVATIVE OF THE SHAPE FUNCTIONS 
C 

DDSHPF ( 1 ) = O.ODO 

DDSHF'F ( 2 > = ( - 6/EL? > T ( 1 2*X/EI .3) 

DDSHPF ( 3 > » DDSHPF (2) 

DDSHPF (4) s O.ODO 

DDSHPF (5) = (-4/ELH (6*X/EL2> 

DDSHPF ( 6 ) = DDSHPF (3) 

DDSHPF (7 ) = O.ODO 

DDSHPF < 8 ) = (6/EL?>-( 12*X/EI. 3) 

DDSHPF (?) = DDSHPF (8) 

DDSHPF (10) = O.ODO 

DDSHPF (11) = (-2/ELH (6*X/EL2) 

DDSHPF (12) = DDSHPF (11) 

C 
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900 CONTINUE 

C 

C 

RETURN 

END 
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SUBROU TINE UFTSEN< PSIE< fNff NELMr lPNAht > 
cp ** ********** *** ******************************************** ** 

CP* 

CP* getsen: branches to the appropriate constraint type. 

CP* 

********* ************ ******* ************************ ******** 

CP* 

c:p* description: 

CP* 

CP* 'GETSEN' BRANCHES TO THE APPROPRIATE CONSTRAIN I TYPE. 

CP* THE AVAILABLE CONSTRAIN I TYPES ARE: 

CP* COMP - COMPLIANCE 

CP* D I SP - DISPLACEMENT 

CP* STRESS 

CP* 

CP************************************************************* 


CP* 


CP* 

PS IB 

CP* 


CP* 

NT 

CP* 

NELM 

CP* 


CP* 


CP* 

IPNAME 

CP* 


CP* 



THE CONSTRAINT VALUE IS KLlURNLD F UR UNITE 

DIFFERENCE EVALUATION 

COUNTER FUR THE FINITE DIFFERENCE 

TOTAL NUMBER OF ELEMENT'S IN FINITE ELFMtNI 

MODEL. (HIS IS IN IFAD DATA BASE r AND I S 

RETRIEVED WHEN DATA BASE IS OPENED. 

NAME OF FINITE ELEMENT ANALYSIS - CHARACTER 
INTERGEK CONVERSION. 


CP* ******** ************************* ******************** * ****** 


c 

INCLUDE 'CAEGSDR. INCJ INPLIC.SPC' 

INCLUDE 'CAEGSDR. INC! ACLIPN.MON' 

INCLUDE ' CAEGSDR . INCH LNlL.MON' 

INCLUDE 7 CAEGSDR* INC3 ELEUES. MON' 

INCLUDE 7 CAEGSDR. INC 3 SVb'CTR.MON' 

C 

COMMON /MEMARY/ IARRY (60000) 

CHARACTER NAME*8» A*4 r B*4 
C 

DIMENSION IPNAME < 2 ) 1 IDBPTK(4> t IERDF V( 2) .PSJ.B'2) 
C 

DATA ISIZE/60000/ 

C 

NT = NT+1 
C 

60 CALL ACCFND ( 7 j IPNAME r 1 r Or Or 1LKR> 

IF( IERR.NE'.O) GO TO 800 
C 


C OPEN THE IFAD DATABASE 
C 


CALL INENTR ( 1 SIZE v IEKBEV? 1 PNAME' > 1 D HP TK r 3 r I S T A T ) 
IF < ISTAT . NL « 0 ) CO TO 810 
A = CHAR ( IPNAME ( I ) > 

B * CHAR ( IPNAME ( 2 ) ) 

NAME * A//B 


C BRANCH TO THE APPROPRIATE CONSTRAINT ROUTINE 
C 

100 IF(ICT.EO.l) CALL COMP ( PS I B f NT > NK I M ) 

I F < I C T . EG . 2 > CALL D I SP (PS IB • N T r NE l M > 

IF < ICT . EG . 3 > CAL L S TRESS < PS IB HKT NLLM r NAME ) 
GO TO 820 
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C WRITE ERROR MESSAGES 10 THE SCREEN 
C 


000 

F’R INI- 

875 f 

IERR 


GO 10 

820 


810 

PRINT 

876 » 

IERR 


GO TO 

820 



c 

c 

C CLOSE IFAD DATABASE 
C 

820 CALL I NEXT. T ( IDBP Tk» I ERR ) 
IF ( I ERR . NE « 0 ) GO TO 830 
C 

GO TO 900 
C 

850 PRIN1 87/ » I ERR 
C 

900 CONTINUE 


i: 

C FORMAT STATEMENTS 


C 

875 

FORMAT ( 1 X r ' ACCFNU 

RETURNED 

WITH 

ERROR 

' tl4) 

876 

FORMAT (IX* ' INEN 1 K 

KE TURNED 

WJ f H 

ERROR 

M4) 

877 

FORMAT (IX? ' INEXIT 

RE TURNE D 

WITH 

ERROR 

' tiA) 

1000 

FORMAT <2A4> 





1001 

FORMAT (F8. 5) 





1 004 

C 

c 

FORMAT (A6) 





c 

RETURN 

END 
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SUBROUTINE LC1 6 < X r Y r Z • XL » YL t ZL ? TB > 

CP*************** ***** * ******** **************** if ****** ******* 

Qp# 

CP* LC16J TRANSFORMS GLOBAL COURU TO LOCAL COORD FOR 'J601' 

CP* 

CP******************************************* ****** ****** ***** 
L P * 

cp* description: 

CP* 

CP* 'LC1 6' TRANSFORMS THE TRIANGULAR (/.TUBAL COORDINATES 

CP* OF ELEMENT 1601 TO LOCAI COORDINATES 

CP* 

CP ********************************** ***********t^t^*.#tt** + *^*^ 


CP* 



CP* 

X 

GLOBAL X COORDINATE 

CP* 

Y 

GLOBAL Y COORDINATE 

CP* 

z 

GLOBAL Z COORDINATE 

CP* 

XL 

LOCAL X COORDINATE 

CP* 

YL 

LOCAL Y COORDINATE 

CP* 

ZL 

LOCAL Z COORDINATE 

CP* 

CP* 

TB 

TRANSFORMATION MATRIX 


CP*******************************************f: <t**t *********** ♦ * 

C 

INCLUDE ' C AEGSDR ♦ INC 3 IMPLIC.SPC' 

INCLUDE ' C AEGSDR . INC. 3 ACCIPN.MON' 

INCLUDE ' C AEGSDR . INC J ELEDES.MON' 

C 

DIMENSION X(3>fY<3> f 2 ( 6 > ? TD < 6 r 6 > ? CDPRL K 3 1 3 ) f>L ( 3 ) fYL<3> > 

* ZL < 3 ) »C00R(3f 3) »V12(3> ? 013(3) *VfR 3) • T ' 3 , 3 > 

(; 

CHARACTER STYPEM 
C 

DATA IREF/l/f ITYPE/1/ 

c; 

C*** INTIALIZE VARIABLES 

C 

C 

DO AO 1*1 r 3 

COOR<lrI> = X<I> 

COOR ( 2 r I ) * y < I ) 

COOR < 3 > 1 ) * Z'.I) 

AO CONTINUE 
C 

C*** GET THE VECTORS PARALLEL TO THE 1-2 AND 1-3 SIDES 

C 

CALL UMVEC(COOR(t r 1 ) * CUUR < 1 ,2) *V12r IERR) 

CALL UMVEC < COOR ( 1 r 1 ) r COOR < 1 » 3 > f VI 3 » IERR ) 

C 

C OBTAIN NORMA! V12 X V13 

c: 

CALL UMVCRS< V12t V13> VNfOf IERR) 

C 

C OBTAIN LOCAL TRANSFORMATION MATRIX I I i 
C 

CALL EUTSCS ( VNf Tf IERR > 

c: 

C OBTAIN THE NODAL TRANSFORMATION MATRIX AS AN ASSEMBLFOCE OF CTJ 
C 

DO 50 K>1?3 
DO 50 J~ 1 * 3 
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TB(JrK) = T(JfK) 

TB< J+3rK + 3> = 1 (J» K ) 

50 CONTINUE 

C 

C GET LOCAL COORDINATES 
C 

CALL UrtXATB ( I » COOR > COORL * 3 1 3 * 3 ) 
DO 60 1*1 » 3 

XL < I ) * COORL (It I) 

YL ( I ) * COORL (2* t > 

ZL ( I ) = COORL < 3 f I ) 

60 CONTINUE 

C 

RETURN 

END 
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SUBROUTINE SKSTQ&<1LCN> 

CP************* ***************************** ****** ************* 

CF'# 

CP* SRST05 S CALCULATES THE ADJOIN I LOAD FOR 1 HF. BK.AH ELEMENl 
CP* 

CP* ************ ******************* ****************t ************ 
CP* 

cp* description: 

CP* 

CP* ' SRST05 •' CALCULATES THE ADJOINT LOAD FOR BENDING 

CP* STRESS IN A 1-DIMENSIONAL BEAM 

CP* ELEMENT AND CREATES A RESTART FILE SO THAT A RESTART 

CP* OF THE FINITE El.Ertt.Nl MOL'EI CAN BE MADE. I HF 

CP* RESULTING DISPLACEMENTS CAN THEN BE USED TO CAI CULATE 

CP* THE STRESS C0NSTRA1NI DESIGN SENSITIVITY. 

CP***#* THIS ROUTINE ONLY WORKS FUR BEAMS LYING IN THE 
CP* X-GLOBAL OR Y-GLODAL PLANE. 

CP* 

GP******************************************* *********** ******* 
CP* 

CP* ILCN INTERNAL LOAD CASE NO. OF ORIGINAL I. OAR 
CP* 

CP******************************** ************* *t* ********* ♦ ft* 

C 

INCLUDE ' CAEGSDR * INC 3 IMPL1C . SPC ' 

INCLUDE ' CAEGSDR. INC 3 ACCIFN.MON' 

INCLUDE ' t AEGSDR* INC 3 CNTL.MON' 

INCLUDE ' CAEGSDR. INC 3 ELUDES * NON ' 

INCLUDE ' CAEGSDR . INC 3 SVECTR.MGN' 

C 

DIMENSION X<2>fY<2)fZ.(2>f DDSHPK ( 1 ? ) f DATN ( 50 > • DUE ( 1 00 ) f 

* SHPF ( 12 >fCFBUK (6) fGPLW<3> fAI GP<12> f CD*. 2 > 6 ) » W lU ( 3 ) f 

* IDIR(2)rAL(12) 

C 

DATA GPLW/-*.774596669241483D0f .ODOr 

* .774596669241483D0/ 

DATA WfW/ . 5555555555555556D0 ? . 088888888888888900 f 

* . 555555555555555 6D0/ 

DATA KT/3/f IREF/1/fMP I /J / 

C 

C INITIALIZE VARIABLES 
C 

ND * ND0F*NUNPE 
DO 10 J=t,ND 
AL ( J ) = 0.0 
10 CONTINUE 
C 

C GET DEPTH AND WIDTH OF THE BE Ah 
C 

CALL ACCEPR < 2 f .1 F’NEPRf IPTAPf Of BUh fLEN f I ERR ) 

IF(IERR.NE.O) GO TO 809 
H = 2*BUF(9) 

B = 2*BUF(10> 

C 

C GET MODULUS OF ELASTICITY 
C 

CALL ACCMAT<2fIPNMAT* NMAT t MP I fHUF f LFNf IERR ) 

IF( IERR.NL.O) v GO TO 808 
E = BUF(b) 

C 

C GET Xf Yf AND Z OF ELEMENT NODES 


c 

CALL ACC EL C ( 2 1 IF'NELC * K INT * IRFF » BUF » I.EN* IERR > 
IF< I ERR* NE •(>> 6(3 TO 907 

M = 1 

DO 200 J=1 f 6 r 3 
K = J+l 
L = J+2 
X < M ) = BUF<J> 

Y (M) = BUF(K) 

Z(M> = BUF (L > 

M = M+l 

200 CONTINUE 

C 

C CALCULATE ELEMENI LENGTH 
C 

PX = X < 2 > - X < 1 ) 

DY = Y(2)-Y(l> 

DZ = Z(2)-Z(l> 

EL = DSQRT ( DX*DX+DY*DY +DZ*PZ ) 

C 

XhP = 1 « DO/EL. 

C 

DO 300 K = 1 r 3 
FSI = GFLU(K) 

CALL EU3DSB<PSI r SHPF t DDSHPF r 2 > EL > 

C 

C CALCULATE ADJOINT LOAD 
C 

ALGP < 3 ) = -.5D0*H*L#XMK*DDSHPF(3> 

ALGP ( 9 ) = -.5D0*H#E»Xf1P#DDSHPF<y) 

IF ( DX . EQ ♦ 0 • AND • DZ * ER « 0> GO TO 250 
ALGP ( 4 ) = O.ODO 

ALGP ( 5 ) = -.5U0*H*E#XMK*nDSHPF<5) 

ALGP<10> = 0 « ODO 

ALGP (11) = -.5H0*H*E*XrtP*DDSHPh <11 ) 

GO TO 255 

250 ALGF* ( 4 ) = -.5D0*H*E*XMk*H1»SHPF<5> 

ALGP < 5) = 0 ♦ ODO 

ALGP (10) = -.5D0*H*E*XrtH*D»SHPF(ll> 

ALGP (11) = 0 * ODO 

C 

C SUM ADJOINT LOAD OVER GAUSS POlNlS 
C 

255 DO 260 J-1»ND 

AL ( J) = AL ( J ) + ALGP(J) 

260 CONTINUE 

300 CONTINUE 
C 

N * -1 

DO 400 J”1 f NUNF'E 
DO 350 K=1>N00F 

IF ( AL ( K + J+N ) « LQ .0*0) GO TO 350 
WRITE (Ilf 864 > IEXTNNC J> fK»AI (K+J+N) 

350 CONTINUE 

N = N + 5 

400 CONTINUE 

C 

C 

GO TO 820 
C 

C WRITE ERROR MESSAGES TO THE SCRFEN 



c 

802 

(307 

(308 

809 

C 

C 

820 

C 

C 

860 

864 

872 

876 

878 

879 
1001 
C 

C 

C 


PRINT 872 f I ERR 
60 TO 820 
PRINT 879 f I ERR 
GO TO 820 
PRINT 879 f TERR 
GO TO 820 
PRINT 876 r I ERR 
GO TO 920 


CONTINUE 


FORMAT ( IX* 'GP~' r I2r 4Xr 'U=>' 9 fr~l 1 « 5» 4Xv 'WXX=' ?E 
*E1 1 ♦ 5» 4X> 'AUXX=' > Ell .b> 


FORMAT (2(I4f2X>fK16*8) 
FORMAT (lXr 'ACCCND RETURNED 
FORMAT (lXf ' ACCEPR RETURNED 
FORMAT < 1 X r ' ACCELC RETURNED 
FORMAT < IX f ' ACCMAT RETURNED 
FORMAT C E12 . 5 ) 


WI TH 

ERROR 

' t }A) 

UJITH 

ERROR 

' tlA) 

WITH 

ERROR- 

' f XA > 

Mil H 

ERROR 

' f 14) 


RETURN 

END 
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SUBROUTINE SRST1 1 ( ILCN f IE ) 

CP*************************************** ************ ********** 
CP* 

CP* SRSTli: CALCULATES THE ADJOIN I LUAD I OR THK PLANE STRESS EL 
CP* 

CP********** ********************************************** ***** 
CP* 

cp* description: 

CP* 

CP* 'SRSTli' CALCULATES THE ADJOINT LOADS FOR A FOUR AND 

CP* EIGHT NODE PLANE STRESS ELEMENT AND THEN CREATES A 

CP* RESTART FILE SO THAT THE FINITE ELEMENT MODEL CAN BE 

CP* RESTARTED SO I HAT THE RESULTING STRESSES AND S ! RAINS 

CP* CAN BE THEN USED 10 CALCULATED THE STRESS CONSTRAINT 

CP* DESIGN SENSITIVITY* 

CP* 

CP************************************************************ 

CP* 

CP* ILCN INTERNAL LOAD CASE NO. OF ORIGINAL LOAD 

CP* IE EXTERNAL ELEMENT NO. CONSTRAINED 

CP* 

CP************************************** *************** ******** 

c 

INCLUDE ' CAEGSDR . INC J IMPLIC.SPC' 

INCLUDE ' C AEGSDR. INC 3 ACCIPN.MON' 

INCLUDE ' CAEGSDR. J NC3 CN1L.M0N' 

INCLUDE ' CAEGSDR. INC!) El EDES ♦ MON ' 

INCLUDE 'CAEGSDR. INC J SVECTR . NON ' 

C 

DIMENSION SHPF < 8) f GPL < 2 f 4 > f SBUK < 50 > f BUT < 3 00 > f 

* DSHRGX ( 8 ) f DSHPG Y ( 3 ) r DSHPL ( 2 f 8 ) f X l 8 ) fY(8)fZ(8)f 

* DG < 3 > i ALGP < 1 6 > f AL < 1 6 > » E < 3 f 3 > f U < 9 ) f C < 3 > f B < 3 f 1 6 > f 

* EMT(l) 

C 

DATA GPL/2*-. 57735027 f . 57735027 f -. 57735027 » 

* 2* ♦ 57735027 f-. 57735027 f .5773502 7/ 

DATA KT/3/ f I REF/1 / f MPT/3 / f I TYPE/1/ 

C 

C** INITIALIZE VARIABLES 
C 

ND = NUN PE* ND OF 
DO 10 J=1f ND 
10 AL ( J) = 0. 

IF(IE.GT.l) GO TO 15 
C 

C** GET X AND Y FOR JACOBIAN EVALUATION 
C 

1 5 CALL ACCELC ( 2 f IPNELC f KIN f f I REF f SUE f LFNB f 1 ERR ) 

IF(IERR.NE.O) GO TO 807 
M = 1 

DO 20 L = 1 f LENBf 3 
J = L + l 
K = L + 2 
X ( M ) * DUFCL) 

Y ( M ) = FUF(J) 

Z ( M ) = BUF <K) 

M = flfl 

20 CONTINUE 

AREA = AREAQ(XfY) 

XMP = 1. DO/ AREA 
C 
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C#* CALCULATE ADJOIN I LOAD 
C 

C*** LOOP OVER GAUSS POINTS 
C 

IC « 0 
H = 1 

DO 80 KK-1 r NSVAL 
J « M+l 
L = M+3 
C 

C*** GET DERIVATIVES OF STRESS FUNCTION G 
C 

C*** VON MISES ♦ (::=<SIGX**2 + SIGY**2~SIGX*SIGY+3*SI6XY**2)**.5 
C 

IF( IC « GT ♦ 0) CO TO 25 

CALL ACCFES(2fIPNFESfK1NTf IREFfII.CNfSBUFfI. LNf I ERR > 
IF ( I ERR t NE ♦ 0 ) GO TO 801 
25 I F < 1ST « EG . 1 ) GO TO 30 

VMS = <SBUF(M)**2 + SBLJM J) **2-SBUM M > *SBUF < J>+3* 

* SBIJF<L>**2)**.S 

DG ( 1 > = • 5*<2*SBIJF<M>-SBUF< J>>/VMS 
DG ( 2 ) = .5*<2*SBUF(J)-SBUF<h)>/VrtS 
DG ( 3 ) = <3#SBUK(L)) /VMS 
GO TO 40 


C 

C *** PRINCIPAL STRESS 
C 

30 THAX = < ( .5*<SBUF<M)-SBUF< J ) ) ) **2+SBUF ( L ) **2 ) ** ♦ 5 

DG ( 1 ) « .5+.25*<SBUK(M>-SBUF< J) )*< 1/TMAX) 

D G ( 2 ) = *5-*25*(SBUF (M)-SBUF ( J) >*( 1/fflAX) 

DG < 3 ) = SBUF < L ) ♦ ( 1/ ThAX ) 

C 

*10 M = M + 4 

IC * ICT1 

PS I = GPL < 1 • KK > 

ETA = GPL ( 2 f KK > 

I F ( I S TYP . t G « 2 ) CALL EU2DLQ < PS J f FTA f K f f SHPF f DSHPL f 

* DSHPGX f DSHPGY f MET Jf X f Y f IERR ) 

IF( 1ST YP. EC). 4) CALL EU2DPG < PSI f E I A f K f f SHPF f DSHPL f 

* DSHPGX f DSHPGY f MET J f X f Y f I ERR ) 

I F ( IERR * NE . 0) GO TO 80? 

t: 

C*** GET ELASTICITY MATRIX E 

C 

CALL EU2DSS ( D * TMP t EM X f NhAT f ITYfiAl 1 1 TYPE ) 

N * 1 

DO 50 JJ=1 f 3 
LL = N+l 
LLL = N+2 
E(JJfl) = D(N> 

E ( JJf 2) = IUL.L ) 

E (JJr 3) = D (LLL > 

N « N+3 

50 CONTINUE 

i; 

C*** GET STRAI N-DISPLACEMEN I MATRIX B 
C 

CALL SD1 1 ( Bf DSHPGX . DSHPGY f NUNPE ) 

C 

C*** CALCULATE £ DG J*CE J*CB3#MP 


C 
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CALL UMXA&atGrE.Cf 1 tSt 
CALL UMXAE< ( C r Br AL. GP > 1 • NlJ > :s ) 

DO 60 J J-l t ND 

60 ALGP ( JJ > = ALGP(JJ>*Xf1P 

C 

C*** SUM ADJOINT LOAD LIVER GAUSS POINTS ( IN TEGRATF UVL R ELEM) 

C 

DO 70 JJ«lrNU 

70 AL(JJ) * AL<JJ>+ALGP<JJ>*DETJ 

80 CONTINUE 

C 

C** WRITE ADJOIN 1 LOAD TO RESTART KILE 
C 

N * -1 

DO 95 JslfNUNPE 
DO 90 K-l 9 NDDF 

WRITE < 11? 2004 ) IPX TNN< J> vKr AL<K + J+N) 

90 CONTINUE 

N = N+l 

95 CONTINUE 

C 

GO TO 820 

C WRITE ERROR MESSAGES 10 THE SCREEN 
C 

801 PRINT 871 » I ERR 

GO TO 820 

007 PRINT 878 t I ERR 

GO TO 820 

009 PRINT 376 ? I ERR 

GO TO 320 

C 
C 

820 CONTINUE 
C 

c 

871 FORMAT < 1 X * ' ACCFES RETURNED WITH ERROR 'f.(4> 

876 FORMAT < 1 X * ' SHAPE FUNCTION ROUTINE RETURNED WITH ERROR '> 
*I4> 

078 FORMAT (IX* 'ACCELC RETURNED WITH ERROR ',M> 

2001 FORMAT ( A ) 

2004 FORMAT ( lXrI4rlX*I2?lX?E16»9) 

C 

c 

RETURN 

END 
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SUBROUTINE SRST16 ( ILCNf IE ) 

CP**-*******#*# ** ***************************** ****************** 

C p* 

CP* SRST16 5 CALCULATES THE ADJOINT LOADS FOR A TRI. BENDING LI.. 

CP* 

CP************************************************************* 

CP* 

cp* description: 

CP* 

CP* ' SRST16 ' CALCULATES THE ADJOINT LOADS FOR A I R1 ANGULAR 

CP* PLATE BENDING ELEMENT AND CREATES A REST ART FILF. SO 

CP* THAT THE FINITE ELFMtNT MODEL. CAN BE RF STARTED* I HE 

CP* RESULTING STRESSES AND STRAINS ARE THEN USED IN THr 

CP* STRESS CONSTRAINT DESIGN SENSITIVITY CALCULATION IN 

CP* THE ' ST1601 ' SUBROUTINE. STRESS TYPES ^PRINCIPAL 

CP* 2 * DON MISES 

CP* 

CP************************************************************* 

CP* 

CP* ILCN INTERNAL. LOAD CASE NO. OF THE ORIGINAL I HAD 

CP* IE EXTERNAL CONSTRAINED ELEMENT NO. 

CP* 

CP************************************************** *********** 

C 

INCLUDE ' CAEGSDR . ) NC J IMPLIC.SPC' 

INCLUDE 'CAEGSDR. INC 3 ACCIPN.MON' 

INCLUDE ' CAEGSDR. INC3 CNIL.MON' 

INCLUDE 'CAEGSDR. INC 3 El PDFS. MON' 

INCLUDE ' CAEGSDR. INC3 SVECTR.MON' 

C 

DIMENSION X<3> fY<3) f Z< 3 > f BUF <50 > fAL<18> • SBUF‘ (SO) ? DO ( 3) t 

* CG < 3 > f XL < 3 > f YL < 3 > f ZL < 3 > f C < 3 > r H < V > f Eh 1 < 27 > f E < 3 f 3 > f 

* TB<6f6> fSIG<3f 3> f ALM( 18) 

C 

C 

DATA IREF/1/fITYPE/1/ 

c 

DO 10 J~ 1 f 1 8 
AL < J ) * O.ODO 
10 CONTINUE 

C 

C** GET Xf Y AND 2 

C 

15 CALL ACCELC ( 2 1 1PMEI CrKINl fIREFf BUFfI ENBfJFKR) 

IF(IERR.NE.O) GO TO 807 
M = 1 

DO 20 L=1fLF.NB»3 
J * L41 
K = L-L2 
X ( M> * BUF<1.) 

Y < M ) - BUF(J) 

Z<M> = BUT ( K ) 

M = M+l 

20 CONTINUE 

C 

C** GET LOCAL COORDINATE SYSTEM 
C 

CALL LC16<X> YfZfXLfYLfZLfI B) 

C 

C* GET ELASTICITY MATRIX CE3 

C 
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CALL EU2DSS < D , I HP ?EHY , NrlAT t 3 TYMAT ? I TYPE > 

N = 1 

DO 25 JJ=1 r 3 
E(JJfl) = D < N ) 

E < J J i 2 ) = D ( N+l ) 

E ( J J f 3 > = D( N+2 ) 

N = N+3 

25 CONTINUE 

C 

C GET THE THICKNESS OF THE PLATE 

C 

CALL ACCEF’R( 2 f lF’NEPR r J P I AP> 0 • PUh *LEN r IERR > 

IF( IERR ♦ NE* 0) 60 TO 909 

THK * BUF ( 25 ) 

c 

C** GET STRESSES AT THE MIDSIDE NODES OF TRIANGLE 
C 

CALL S S M D 1 6 < X L * Y L f I L C N t S I G ) 

C 

C** CALCULATE THE ADJOINT LOAD WHfcRE 0^CDG3 * CE J*C BJfMr* T/2 
C** AL (18) - < TARE A/ 3 ) # < Q ( 0 f *5? *5 )+Q < ♦ 5 f 0 r ♦ 5 ) + CM ♦ 5 f « 5 f 0 ) 

C 

TAREA = EUTRI A ( X ? Y ) 

DO 80 I T= 1 ? 3 
C 

C*** CALCULATE CDGT THE DERIVATIVE VECTOR 
C 

IF ( 1ST « EG« 1 ) GO TO 30 
C 

C *** VON MISES: G=SGRT ( SI GX#*2-f SIGY**2-S1GX#S J(.'Y*M*SIGX 2 > 

C 

VMS = DSGRT<SIG<ITf1)**2+SIG< IT f 2>**2-SIP.< ITf 1 ) *SIG< IT* 2) +3* 
* SIG<ITf3>**2> 

DG ( 1 ) - .5D0*<2.D0*SIG< IT, I )-SIG(ITf 2) )/Vf1S 
DG ( 2 ) = .5D0*<2.D0 *SIG<ITf2>-SIG<ITf J >>/VMS 
DG ( 3 ) = <3.D0*SI6< IT, 3) )/VMS 
GO TO AO 
C 

C** PRINCIPAL STRESS 
C 

30 THAX = DSQRT < < *5D0# < SIG < ITf 1)-*SIG<I T f 2 ) ) ) **2 f SI C ( 3 T f 3 ) **2 > 

DG < 1 ) * « 5D0+ ♦ 25D0* ( SI G ( I T * 1 ) -S 1 1» O T f 2 ) > / I MAX 
DG < 2 > = • 5D0- . 25D0* < SIG ( IT f 1 > -SIG ( I T * 2 ) 3/1 MAX 
DG < 3 ) = SIG< IT r 3 ) / fMAX 
C 

C*** CALCULATE 1X3 * CDG3 * CE3 
C 

AO CALL UMXAB(UGfE»Cf 1 f3f3> 

C 

C*** CALCULATE CAL3 * CC3 * CEO * T/2 t XMP 
C 

CALL AL1 6 < XL * Tl. » C t ALM t I Bp THK r IT) 

DO 50 K— 1 f 18 

AL ( K ) = AL ( K ) + (TAREA/3 *ODO)#ALM(K) 

50 CONTINUE 

BO CONTINUE 

C** WRITE ADJOIN f LOAD TO RESTART FILE 
C 

N = -1 

DO 95 J* 1 9 HUN PE 


DO 90 KMrNDOF 

IF (ALCK+J+N) .FG.0.01'0) (>0 TO 90 

WRITE < 11 >2004 ) 1EXTNN < J ) f K f AL <K+ J IN) 
90 CONTINUE 

N ■ N+5 

95 CONTINUE 

C 

GO TO 820 

C WRITE ERROR MESSAGES TO THE SCREEN 
C 

801 PRINT 871 » I ERR 

GO TO 820 

807 PRINT 878 r IERR 

GO TO 820 

809 PR IN I 879 ? IERR 

GO TO 820 
C 
C 

820 CONTINUE 

C 

c 

871 FORMAT ( IX f ' ACC FES RETURNED WITH ERROR * • J 

878 FORMAT (IX* 'ACC EEC RETURNED WITH ERROR ' f ) 

879 FORMAT (IX*' ACCEF’K RETURNED WITH ERROR 'rl 

2001 FORMAT ( A > 

2002 FORMAT ( IX r 14 f2Xf 'DISP- 7 > El 6, V) 

2004 FORMAT <1XfI4f1XfI2f1X»E16.9) 


RETURN 
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SUBROUTINE SShJLU6<X r Y* 3 LCNf S llO 

CP*********** ****************** ******************************** 

Cp* 

CP* SSMD1 6 J CALCULATES STRESSES AT THE MJDSJDF HOPES OF A TKI . 

CP* 

CFv* ***** *********************************************** ******** 

cp* 

cp* description: 
cp* 

CP* 'SSMD16 ' CALCULATES STRESSES AT THE MIDSIDE NODES OF 

CP* TRIANGLE 1601. CSIG3 « LE3* C BD # C D IS 3* T/2 IS t Hi- 

CP* FINITE ELFMENI FORMULA FOR CALCULATING THE STRESS. 

CP* THE STRAIN-DISPLACEHH N I MATRIX CB3 IS LOCATED )H 

CP* COLUMNS Ay S? AND 6 OF THE CW3 MATRIX . LET IS THE 

CP* ELASTICITY MATRIX* 3X3 FOR ISOTROPIC MATERIAL* AND 

CP* CDIS3 IS THE DISPLACEMENT VECTOR* WHERE THE IOSPLS. 

CP* ARE TAKEN AT THt NODES OF THE TRIANGLE. BECAUSE 

CP* IFAD CALCULATES THE STRESS RESULTANiS* THt T/2 TERM 

CP* MUST BE INCLUDED IN I HIS CALCULATION. 

CP* 

C P************ ************************************ ************* 
t P* 

CP* X THE LOCAL ELFMENI* X COORDINATE 

CP* Y THE LOCAL ELEMENl Y COORDINATE 

CP* ILCN THE INTERNAL I OAD CASE NUMBER FOR THE ORIGINAL 

CP* ORIGINAL MODEL 

CP* SIG THE MIDSIDE NODE STRESSES - RETURNED VALUE 

CP* 

CP* ****** **************************************** ************** 

c 

INCLUDE ' C AEGSDR . INC 3 IMPLIC.SPC' 

INCLUDE ' C AEGSDR. INC 3 ACCIPN.MON' 

INCLUDE 'C AEGSDR. 1NCJ CNFL.hON' 

INCLUDE M. AEGSDR. INC 3 EL.EDES . MON ' 

c: 

EQUIVALENCE <NDAT<14> fIPK) 

0 

DIMENSION X(3> r Y(3> fGFTSC3f3) rU ( 18 , 7 ) , CFhVt </> ?VJK3,/> r 

* DIS(9>fD(9 ) fEC3f3> fEMT <27) rPUM 100) rLtf( 3 f9> » 

* SIG(3f3>fXL<6>fYL<6> 

C 

DATA GPTS/O.ODOr .5D0 f .5D0» .bDOrO.ODOr .M.iOf .5D0f .SDOrO.ODO/ 
DATA I TYPE/1/ 

C 

C** GET DISPLACEMENTS 
C 

DO 10 JbIfNUNFE 

CALL ACCCN)K2r IFNCNDr INTNN( J) rl » )LCNfCFBUF fL.ENf IERR) 
IF(IERR.NE.O) GO TO 802 
DO 10 K*lvLEN 

C D < J f K > * CFBUMK) 

.1.0 CONTINUE 

M « 1 

DO 20 J=1fNUNPE 
DIS< M> » CD( Jf 3) 

DISC M+l ) » CD< Jf 4 > 

DISC M + 2 ) = CDCJfS) 

M « M f 3 

20 CONTINUE 

C 

C** GET ELASTICITY MATRIX 
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CALL EU2DSS ( D . TMR .EH ft NMAT . 1 I YhAT f IT YPE ) 

N = 1 

DO 30 JJ= 1 ? 3 
E ( J J * 1 ) * D(N) 

E ( JJf 2> » D(N+1> 

E ( J Jr 3 > * D<N+2> 

N = N+3 

30 CONTINUE 
C 

C** GET MATERIAI THICKNESS 

C 

CALL ACCEPR< 2. 1PNLPR » IP TAB * 0 . DUF . LENr I ERR) 
IF(IERR'NE.O) GO TO 80? 

THK * PUP (25) 

C 

C ** GET LOCAL X AND Y COORDINATES 

C 

CALL MOVESP<XLf X.3*IPR> 

CALL M0VESP<YL»Yr3*IPR> 

C 

C** CALCULATE STRESSES AT THE MIDSIDE NODES 
C 

DO 100 I NT =1*3 

C 

C*** GET SHARE FUNCTION AT MIDSIDE NODE 
C 

CALL SF 1 501 < XL . YL . GR I S ( 1 . ] N f ) .W. ED) 


C 

C*** CSIGJ = CED*CD3*CDIS3*T/2 


C 

DO 50 M“1 . ? 

DO 50 K = 1 . 3 
GASH = O.ODO 
DO *10 J- 1 * 3 

GASH = GASH + fc ( K . J ) * W ( M . J>3 ) 

40 CONTINUE 

EB(KfM) = GASH 
50 CONTINUE 

DO 60 II = 1.3 

SIG ( I N I fll) = O.ODO 
DO 60 M2 = 1.? 

S I G ( I N f . 1 1 ) = SIG(INT.I1)+EB(I1. M2 ) *D1S < M2 > 
60 CONTINUE 

DO 70 12=1.3 

SIG( INT.I2) = SIG(INf.I2) #THK/2 . ODO 
70 CONTINUE 

100 CONTINUE 
GO TO 820 


C 

C 

C** WRITE ERROR MESSAGES TO THE SCREEN 


C 

802 


80? 


C 

C 

820 


PRINT 

872. 

I ERR 

GO TO 

820 


PRINT 

87?. 

I ERR 

GO TO 

820 


CONTI 1 

NUE 
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c 

FORMATUXr '"ACCCNU RETURNED- WITH ERROR 'r)4> 
FORMAT (lXr ' ACCEFR RETURNED WITH ERROR 'tXA> 
i: 

C 


RETURN 

END 



SUBROUTINE SNrtiiil 6 ( X r Y » ILCN r EF'N ) 

CP*$t ************** ****************************** ************** 
CP* 

C P * SNMD16 1 CALCULATES STRAINS AT 1 HE MIDSIDE NODES OF A 'I K I • 
CP* 

CP***************#**#*********^#******#**#***************#***** 

CP* 

cp* description: 

CP* 

CP* 'SNHD16' CALCULATES STRAINSES AT THE MIDSIDE NODES OF 

CP* TRIANGLE 1601. THE FORMULAS FOR STRAIN ART J 


CP* 

CP* EX = J/E*<SIGX-V*SIGY> 

CP* EY = 1/E* < SIGY-V#S1 GX ) 

cp* exy * aiuxr/o 

cp* 


cp************ ******************************************* ****** 


CP* 

CP* X 
CP* Y 
CP* ILCN 
CP* 

CP* EPN 
CP* 


THE LOCAL ELEMKN I X COORDINATE 

THE LOCAL ElEMtNl V COORDINATE 

THE INTERNAL LOAD CASE NUMBER FOR THt ADJOINT 

LOAD 

THE MIDSIDE NODE STRAINS - RETURNED VALUE 


CP************************************************************* 


INCLUDE ' CAEGSDR * JNC 3 IMPL1C.SPC' 
INCLUDE ‘ CAEGSDR ♦ INC 3 ACCIPN.MON' 
INCLUDE 7 CAEGSDR ♦ INCH CNTL.MON' 
INCLUDE ' CAEGSDR « INC J ELEDES * MON' 


C 

DIMENSION X(3)fY<3) fGnS(JrJ) rU( 18 r 7) »CFBUK<7) rCD<3>/)> 
* DISC 9 > rEMI (27) rBUF(lOO) »EPN(ir3> >SIG(3 t 3) 

C 

DATA MPT/1/ 


C 

Ct* GET CONSTANTS 
C 

CALL ACCMAT ( 2 > I PNMAT t NMAT t MP1 r HUH t l.EN» I ERR ) 
I P < I ERR • NE * 0 ) GO TO 808 
E = BUMS) 


V = DUE (7) 

G * E/ ( 2 • DO* ( 1 * DOT V > ) 


C** GET STRESSES AT THE MIDSIDE NODES 
C 

CALL SSMD16 ( X > Y ? ILCN* S JO 
C 

C** CALCULATE STRAINS AT THt MIDSIDE NODES 
C 

DO 100 IT*lr3 

EPN < I T r 1 ) « 1/E*<SIG<IT. J)-V*SIG<)Tf2>> 
EPN ( I T f 2 ) = 1/E* < S) G< IT f 2) -V#SIG< IT > 1 ) ) 
EPN(IT»3) » SIG(ITr3>/G 
100 CONTINUE 

GO 10 820 
C 
c 

C** WRITE ERROR MESSAGES TO THE SCREEN 
C 
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808 PRINT 878, IERRR 
GO TO 820 

20 CONTINUE 

78 FORMA K IX, 'ACCMAT RETURNED WITH ERROR ',14) 

RETURN 

END 



122 


SUBROU r I NE S TRES05 < N T > I r L * PSI BB > 

CP* ************************************** ************** ******** 
CP* 

CP* STRES05J CALC 4 STRESS CONSTRAINT AND SENSIT. H OR A BF AN 
CP* 

CP* ******** ********************************** ************ ****** 
CP* 

cp* description: 

CP* 

CP* ' STRES05 ' CALCULATES THE STRESS CONSTRAIN! AND THE 

CP* DESIGN SENSITIVITY FOR THE 1-D BEAM ELFMLNl IN 

CP* BENDING ♦ INCLUDES TORSION* 

CP* 

CP* *************************************** ********************* 
CP* 

CP* Nr COUNTER FOR FINITE DIFFERENCE 

CP* I EXTERNAL ELFMFNl' NO. FOR ELL MEN I* BEING PROCESSED 

CP* L LOAD CASES - EXTERNAL 

CP* PS IBB STRESS CONSTRAIN! - RETURNED VAI UL“ 

CP* 

CF’ ******************************************* ****************** 

C 

INCLUDE ' CAEGSDR . INC3 IMPL1C.SPC' 

INCLUDE * CAEGSDR * I NCI ACCIPM.MON' 

INCLUDE ' CAEGSDR ♦ J NCI CN f’L . MON ' 

INCLUDE ' CAEGSDR. INC3 EUvDLS.MON' 

INCLUDE ' CAEGSDR . INC 3 SVLCTR . MON ' 

COMMON/L CSDES/DLCSC 90) 

C 

EQUIVALENCE (NDATC14) fIP) f <NDAT<98> f 3 DBL > 

C 

DIMENSION X(3) rY<3) f 7 < 3 ) f DDSHPF < 1 2 > f DATN ( 50 > f!.(2) fBUk (100) f 

* SHPF(12> fCFBUF< 6 ) fGPLU(3) fAI D ( 7t 6 ) f CD ( 2 f 6 ) f W TW ( 3 ) » 

* PS IBB (500) fC ( 6f 2) fD(6f 2) ?T(3f 3) f TB ( 6 f 6 ) f CDL < 2 f 6 ) f 

* ALDL ( 2 f 6 ) f CUOF* ( 3 ? 3 ) 

C 

DATA GPLU/-. 7745966691-MI A 183D0f .ODOf 

* . 7745966692 41 4 J 83D0/ 

DATA WTU/ .5555555555555556D0f ♦ 8888888888888889D0 f 

* . 5555555555555556B0/ 

DATA KT/3/f IREF/l/r MPT/1/ 

C 

DPS I BG = O.ODO 
DPS I HG - O.ODO 
PSIBB(I) = O.ODO 
IF(I.GT.l) GO TO 50 
C 

C GET AREA MOMENT OF INERTIA ABOUT Y-AXIS 
C 

50 CALL ACCEPR< 2r IPNEPR f IPTABf Of BUh fLENf IERR ) 

IF(IERR.NE.O) CO TO 809 
Yl « BUK ( 5 ) 

H * 2*BUF < 9 ) 

B * 2*BUF ( 10 ) 

BH(ND = H 
BW( NT > = B 
C 

C GET HEIGHT DENSITY AND MODULUS UF ELASTICITY 
C 

CALL ACCMAT ( 2 f I PNMAT f NMAT f MPT f BUh f I EN f IERR ) 

IF(IERR.NL'.O) GO TO 808 
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GAMMA * 0*ODO 
E * 30500000 ♦ 170 
V « * 3 DO 

G = E / < 2 ♦ DO * (1 • DO + V ) ) 

C 

C GET INTERNAL LOAD CASE NUMBER FOR ORIGINAL LOAD 
C 

60 CALL ACCLCSC It IFNLCS p IDBLp 1 t 0 1 1EF<R ) 

IF < I ERR * NE « 0 ) CO TO 805 

CALL ACCLCS(2t lPNI.OSf L(1 ) t 2r DLCSf IERR) 

IF ( IERR » NE » 0 ) GO TO 805 
IL1 = DLCS( 21) 

C 

C GET DISPLACEMENTS AT ElEMhNT ENDS 
C 

DO 70 J*1pNUNPE 

CALL ACC C Nil V/.pI PNCND * INTNN < J ) j 1 f IL 1 p CFBUF p I EN > I ERR ) 
IF ( IERR « NE • 0) GO I 0 802 
DO 70 K=1 ? NDOF 

CD(JpK) = CFBUF <K> 

70 CONTINUE 

C 

C BYPASS ADJOINT LOAD CALCULATION 
C 

IF<NT.GT.1> GO TO 150 
C 

C GET INTERNAL LOAD CASE NUMBER FDR ADJOINT LOAD 
C 

CALL ACCLCS( 1 p JPNLCSp IDBLp 1 pOp IERR) 

IF< IERR . NE « 0) GO TO 805 

CALL ACCLCS(2pIPNLCSpL(2) p2pHLCSp IERR) 

IF ( IERR ♦ NE . 0 ) GO TO 805 
ILCN = DLCS<21) 

C 

C GET DISPLACEMENTS AT ELEMENT ENDS FOR ADJOINT LOAD 
C 

CALL ACCCNTKIp IPNCNUp IDBLp 1 p ILCNpOpOp IERR) 

IF < IERR* • NE • 0 ) GO TO 802 
DO 100 J=1 f NUNPt 

CALL ACCCNLK2p IPNCNUp IN TNN ( J) * 1 p ILCNpCFBUi- pL.ENp IERR) 
IF(IERR.NE.O) 00 TO 302 
DO 100 N= Ip NDOF 

ALD(JpK) * CFBUF<K) 

100 CONTINUE 

C 

C*********** *************************************** ****** <** 

C EVALUATE DISF’LS, AND CURVATURE AT THE GAUSS POINT USING 
C SHAPE FUNCTIONS - ONE PT. FOR CUKV. f THREE PT« FOR BISPL 
C************************* ********************************* 

C GET X> Y f AND Z OF ELEMENT NODES 
C 

150 CALL ACCELC(2pIPNELCpKINTpIREFpBUF p|. ENpIERR) 

IF C IERR. NE ♦ 0) GO TO 807 
M - 1 

DO 200 J=1p?p3 
K * JF1 
LL » JF2 
X ( M ) = BUF(J) 

Y(M) » BUF(K) 

Z(M) = BUF(LL) 
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M = M + l 

200 CONTINUE 

DO 210 J=1f3 

COOR(IfJ) = X(J) 

COOR ( 2 f J ) = Y < J ) 

COOR ( 3 f J> = 2 ( J ) 

210 CONTINUE 

C 

C FORM THE ELEMENT LOCAL COORDINATE SYSTEM FOR DISPLACE MEN t S 
C 

IN3 = INTNN < 3) 

( CALL EUBTM<IN3fBETAfCOORfT r IFRR) 

CALL ZEROSPC IBf36*IP> 

DO 220 J=1 • 3 
DO 220 K”1 f 3 

TB(JfK) = T(JfK) 

TB<J+3fM3> = TB(JfK) 

220 CONTINUE 

CALL UMXABT(iBfC0fCf6f2f6) 

CALL UMXABT < TB r ALD r Dr 6 r 2 » 6 > 

DO 230 J=1f2 
DO 230 K=1 f 6 

CDL. ( J » K) = (.: ( K f J ) 

ALDL(JfK) = IKKfJ) 

230 CONTINUE 

C 

C CALCULATE ELEMENT LF_NGTH 

C 

DX = X ( 2 ) -X< 1 ) 

DY = Y (2 ) - Y ( 1 ) 

DZ = Z(2>-Z(l) 

EL = DSQRT<DX*DX+DY*DY+DZ*DZ> 

C 

C CHANGE LOCAL Y-ROTATIUN FROM POSH IVE TO NEGATIVE IF 
t: BEAM LIES ALONG THE X GLOBAL AXIS 

C 

IF<DX. LT. 0,001. AND. DX.GT. -0.001) GO TO 245 
DO 240 J=1fNUNPE 

CDL ( J f 4 ) = -CDL ( J t 4 ) 

CDL ( J f 5 > = -CDL. ( Jf 5 ) 

ALDL ( J f 3 ) = “ALDL(Jf3) 

240 CONTINUE 

C 

C CALCULATE THE TWISTING ANGLES 

C 

245 UXY = DABS ( ( CDL < 2 f 4 ) -CDL < l f 4 ) ) /EL ) 

AWXY = DABS ( < AL DL ( 2 • 4 > -ALDL < 1 f 4 ) ) /EL ) 

C 

C EVALUATE SHAPE FUNCTIONS FOR DJSPL. - THREE POINT (iUADKATUKE 
C 

B2 = B#B 
B3 * B2*B 
B4 ~ B2*B2 
H2 = H*H 
H3 = H2*H 
C 

DO 300 K=1 r 3 
PSI = GPLWCK) 

CALL EU3DSB ( PS) f SHPF f DDSHPFf 2f EL) 

W = <SHPF<3)*CDL<lr3)+SHKF<5)*CDL(if5)+SHPF<9>* 

* CDL < 2 f 3 ) +SHPF (11) *CDL < 2 f 5 ) ) 
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AW = <SHPF(3)*Al.DI. < J . 3 > 4SHPF ( 5 > * Al DL < 1 . 5 > 4 SHPF < V > * 

* AL.DL <2 . 3) 4SHPF Ci 1 ) #ALDL < 2? 5 ) ) 

C 

C EVALUATE SHAPE FUNCTIONS FOR CURV. - THKEF. FOJNl QUADRATURE 
C 

WXX = <DDSHPF<3)*CDL<if 3 > 4HUSHPF ( 5 ) *CDL < 1 .5)4 

* DDSHPF < 9 > *CBL < 2 f 3 > 4BHSHPF (1:0 *CDL < 2 . 5 > ) 

AWXX = ( DDSHPF ( 3 > *ALDL ( 1 . 3 ) 4)iDSHPF ( 5 ) * AI.DL <1.5)4 

* DDSHPF ( 9 ) * ALDl. (2.3) 4 DDSHPF (11) *ALDL ( 2 . 5 ) > 

C 

C CALCULATE SENSITIVITY VECTORS 
iZ 

XMP * 1 ♦ DO/EL 

IF< NT » 6 T • 1 ) GO TO 250 

STERh = « 5D0*XMP* E* WXX 

PJB « H3/3 # DO- • 42W0#B# < H24B4/ (4«D0*H2>) 

PJH = B*H2-.42D0*B2*(H-B4/< 12. DO*HJ> ) 

DPSIB6 = DPSIBG4-<-GAMHA*H*AU-< fc*H3/12.M0)*AWXX*WXX- 

* PJB*G*WX Y* AWX Y)*UTW(K)*(LL/2« DO ) 

IF (I .NE.ICE(NC) ) GO TO 2 47 

DPSIHG = DPSIHG4<-GAMMA#B*AW-<3.D0*E*B¥H2/12*D0>*AWXX*UXX 

* -STERM~PJH*6*WXr*AWXr>*WrW<K>*(EL/2.D0> 

GO TO 250 

247 DPSIHG * DPSIHG 4 ( -GAMhA*B*AW~ ( 3 . DO*E*B*H?/ 1 2 . DO ) * AWXX 

* *WXX-PJH*G*UXr*AWXY>*WrW(K)*<FL/2.D0> 

C 

250 IF < I ♦ NE ♦ ICE (NO ) GO TO 300 

C 

C CALCULATE PSI<B> - INTEGRAL UP STRESS FUNCTION TOR ELEhfc N f 
C 

STRESS * -*5B0*H#E*WXX*XMP 

PSIBB(I> * PSIBB(I) 4 STRESS*WTW < K ) # ( EL/2 ♦ DO ) 

300 CONTINUE 

I F ( N T « G T ♦ 1 ) GO TO 820 
DPSIB(I) = DPS IBS 
DPSIH(Z) * DPSIHG 
C 

GO TO 320 
C 

C WRITE ERROR MESSAGES TO THE SCREEN 
C 

802 PRINT 872 . I ERR 

GO TO 820 

805 PRINT 875. I ERR 

GO TO 820 

307 PRINT 878. I ERR 

GO TO 820 

(308 PRINT 879. I ERR 

GO TO 820 

309 PRINT 376. I ERR 

GO TO 820 

C 

C 

820 CONTINUE 

C 

c 

850 FORMAT < IX » '***ADJOINl' LOAD IS APPLIED AT EIEMENT'.H) 

851 FORMAT (/.IX.' BEAM WIDTH B= ' . F8 . 5 . 2X > ' BEAM DEPTH® ' . F8 . 5 . 2X 
*, 'E=' .E9* 3.2X. 'IYY=' . E9 . 3 . 2X . ' G AMMA= * .F6.5. ' APPLIED FORCE 
*=' »F8.5> 

(355 FORMAT ( IX. 'NODE = ' . 1 2 . 2 X . ' X= ' . F 1 2 ♦ 5 . 2X . ' Y® ' . E 1 2 ♦ 5 . 2X » ' /= ' f 
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*E12.5,2X, 'RX=' ,E12.5,2X, 'KY=' ,E12.5, 2X, 'RZ=' ,E12.5) 

860 FORMAT < IX, ' GR= ' , 1 2 , 4X , ' W= ' » El 1 . 5 , 4X » ' WXX* ' ,E1 1 . 5 , 4X , 'AW=' 
*E11 .5»4X, ' AklXX* ' » til .5) 

( 36 4 FORMAT < IX, ' NODE* ' , ) 2 , 2X , ' AX= ' »E12.5,2X» ' AY* ' ,H 2. S»?X» 

* ' AZ= ' , E12 . 5 , 2X , ' ARX* ' , E12 . 5 » 2X, ' ARY* ' , E12 . 5 , 7X , ' ARZ= ' , 
♦E12.5) 

872 FORMAT ( IX, 'ACCCND RETURNED WITH ERROR ',I4> 

375 FORMAT < IX* 'ACLLCS RETURNED WITH ERROR ',14) 

876 FORMAT ( IX, 'AUCEPR RETURNED WITH ERROR ',14) 

978 FORMAT ( IX, 'ACCELC RETURNED WITH ERROR ',14) 

37? FORMAT < IX, 'AUCHAT RETURNED WITH ERROR ',14) 

1001 FORMAT < El 2. 5) 

1004 FORMAT (14) 
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SUBROUTINE STRESi 1 ( NT tlfLt F’SI BT ) 

CF : ****************************************** ****** ************** 
CP* 

CF** STRESI IS CALC. STRESS CONSTRAINT AND SENSIT. - PLANE STKfcSS 
CP* 

cp*********************** ****** *********************** ********* 
CP* 

CF** description: 

CP* 

CP* ' STRESI 1 ' CALCULATES THE STRESS CONSTRAINT AND IKK 

CP* DESIGN SENSITIVITY FOR A FOUR OR AND EIGHT NODE PLANE 

CP* STRESS E! ILMEN f WITH TRACT I ON ? SELF WEIGHT NOT INCLUDED 

CP* 

CP* ***************************************************** ******* 
CP* 

CP* NT COUNTER FOR FINITE DIFFERENCE 

CP* I EXTERNAL ELEMENT NO. BEING PROCESSED 

CP* L EXTERNAL LOAD CASE NOS. 

CP* PSIBT . STRESS CONSTRAINT 

CP* 

CF‘**** ********************* ************************************ 

C 

INCLUDE ' C AEGSDR. INCH IMPL1C . SRC ' 

INCLUDE ' C AEGSDR. IHC3 ACCIPN.MON' 

INCLUDE ' t AEGSDR. INC J CNTL.MON' 

INCLUDE 'C AEGSDR. INC: ELEDES. MON' 

INCLUDE ' C AEGSDR. INC 3 SVECTRthON' 

COMMON /LCSDES/DLCS ( 90 ) 

C 

EOUI VALENCE ( NDAT ( 98 > * IDBL ) 

D I MENS I ON X(8)»Y(8)fZ(8>» SHPF ( H ) f GPL < 2 > 4 ) p L < 2 > > ILCN< 2 ) t 

* DATN(SO) *EUF< 100) p PSIBT (500) p SBUF < 50 ) r CF BUK <6) p 

* BF( 4 p 48) i EM I ( X) , DSHPGX ( 8 ) p EGBUF ( 50 ) , UHHPCY (H), 

* DSHF’L ( 2f 8 ) p SIGMA < 6 p 4 > p EPSLN ( 6 p 4 ) p B ( 3 p .i 6 ) p SE < 500 ) • 

* DGt3) pALGPC 16) p AL < 16 ) ? E ( 3 p 3 ) r D ( 9 ) * C ( 3 ) p FS I BPE < 500 > 
CHARACTER YESNOU 

C 

DATA GPL/2*- .57735027 p . 57735027 p -• . S//3t»027 » 

* 2*. 577350:*7r-. 57735027 p .57735027/ 

DATA KT/3/ f I REE/ 1 / p MP1 /> /? .1 TYPE/ 1 / 

C 

C 

SE ( I ) = 0.0 
PSIDPE(I) =0.0 
IF(I.NE.l) GO TO 100 
10 DO 50 JJ=1 r 2 

c: 

C GET INTERNAL LOAD CASE NUMBER 
C 

120 CALL ACCLCSdpjPNI CSp ) M*L p 1 fOp IERR) 

IF ( I ERR . NE ♦ 0 ) GO TO 805 

CALL ACCLCS(2pIPNLCSpL< JJ) p2p DLLSp 1ERR) 

IFUERR.Nf.O) CO TO 805 
ILCN < J J > ■ DLCSC21) 

50 CONTINUE 

C 

C SETUP POINTER FOR STRESS-STRAIN BUFFER 
C 

100 DO 165 

CALL ACCFESd pIPNFESp IDBLp 1 p ILCN< JJ) f Of Op IERR) 
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IF < I ERR • NE * 0 ) GO TO 801 

C 

C GET ELEMENT STRESSES AND STRAINS 
C 

CALL ACCFES<2> IPNKt.S? KINT » XRbK» ILCN( JJ> ? KBUF » LEN ? I ERR > 
I F < I ERR * NE ♦ 0 ) GOTO HOI 
M * 1 

IF (JJ.EQ.2) GO TO 140 
DO 130 K“1 f NSVAL 
J = M-41 
LL “ M+3 

SIGMA ( 1 ? K ) = SDUF(H) 

SIGMA ( 2 f K > = SttUMJ) 

SIGMA (3fK> * SBUF ( I..L > 

M - M*M 

130 CONTINUE 

GO TO 160 
1 40 h * 17 

DO 150 K- 1 f NSVAL 
J = M+l 
LL * M+3 

EPSLN < 1 * K ) = SBUK(M> 

EPSLN<2rK> = SDUF(J) 

EPSLN ( 3 1 K ) = SBUF < LL ) 

M « M + 4 

150 CONTINUE 

160 IF < JJ ♦ EG * 2 > GO TO 170 

IFCNT.GT.l) GO TO 170 
165 CONTINUE 

C 

C GET X AND Y FOR JACOE-I AN H.VAI. NATION 
C 

170 CALL ACCELC < 2 1 IPNEL.C fKINT 1 1REF * BUT fI.ENXf JFRR) 

IF < I ERR « NL * 0 ) GO TO 807 
M = 1 

DO 200 J= 1 1 LEND ? 3 
K *= J+l 
LI. = J+2 
X < M ) * BUF(J) 

Y (M) * BUKOO 
Z ( M ) = BUFCLL ) 

M = h+1 

200 CONTINUE 

AREA = AREAQ(XfY) 

XMP = 1 . DO/ AREA 

c; 

C CALCULATE FORCES AT THE GAUSS POINTS 
C 

DO 250 J = - 1 t NDOF 
DO 250 K~1 1 NSUAI . 

BF(KfJ) * 0.0 
250 CONTINUE 

C 

C LOOP OVER THE GAUSS POINTS 

C 

DO 300 K«1 f NSVAL 
PS I * CF’L ( 1 f K > 

ETA « GPL < 2 r K ) 

C 

C EVALUATE SHAPE FUNCTIONS AT THE GAUSS P01NIS 
C 
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IF( IS TYP ♦ EG ♦ 2) CAI L EU2DL.Q < PSI f ETA f K I f SHPF r NSHK. » 

* DSHPGX f DSHPGY f DET J f X f Y f XERk * ) 

IF< ISTYP.EQ.4) CALL EU2DPG ( PSI fETAfK iVSHkE fDSHkL f 
t DSHPGX f DSHPGY f DET J f X f Y f i FKK ) 

IF( IERR.N£.0> GOTO 809 
300 CONTINUE 

WRI TE < 1 0 f 855 ) 

DO 320 K=1fNSVAL 

320 WK'I TE < 1 Of 854 > Kf < SIGMA < J f K > f J= 1 f NSIG > 

IF < NT . GT 4 1 ) GO TO 34b 
WRI T'E ( 1 0 f 360 ) 

DO 330 K-lrNSVAL 

330 WRITE ( 10r 854 ) Kf <EPSLN< JfK) t J~1 f NSIG) 

DO 340 J-1fNSIU 
DO 340 K=1fNSVAL 

SE ( I ) » SF ( I ) + SIGMA ( JfK)*EPSLN< JfK>*UEU 
*340 CONTINUE 

C 

C CALCULATE SENSITIVITY VECTOR 
C 

DPS IT < I > = - SE(I) 

C 

345 IF ( I • NE • ICE ( NC ) > GO TO 820 

C 

C* CALCULATE PSI<B> - INTEGRAL OF STRESS FUNCTION G V OR ELFhENT 
C 

I F < I ST ♦ EQ ♦ 1 ) GO TO 360 
HO 350 K=1fNSVAL 

VMS = ( SIGMA ( 1 1 K ) ♦♦2+SIGMA < 2 f K ) ♦♦2-SIGMA < 1 f K ) * 

* SIGMA ( 2 f K > +3*SIGMA < 3 f K ) **2 > * * ♦ 5 

350 PSIDPE(I) * PS1PPE ( I > + VMS*DLTJ 

GO TO 300 

360 DO 370 K=lrNSVAL 

TMAX = ( ( * 5*(SIGMA( I f K) -SIGMA (2f K ) ) >**2+SI6MA ( 3f K 
t %% * 5 

370 PSIBPE < I > * PSIBPE< I > + ( .5*<S.|GNA< 1 f K > +SI GMA < 2 f K > )+TMAX> 

* *HETJ 

380 PS IDT ( I ) = PSIDPE < I ) #XnP 

C 

C 

GO TO 820 

C WRITE ERROR MESSAGES TO THE SCREFN 
C 


800 

PRINT- 

870 f 

IERR 


GO TO 

820 


1:101 

PRINT 

871 - 

IERR 


GO ro 

820 


805 

PR INI 

875 f 

IERR 


GO TO 

820 


807 

PRINT 

878 f 

IERR 


GO TO 

820 

- 

80? 

PRINT 

876 f 

IERR 


GO TO 

820 



C 

C 

820 CONTINUE 

c: 
c 

850 

851 
854 


FORMAT ( IXf '*** ADJOIN I LOAD IS APPLIED A T Kl. EMLN I ' f 14 ) 
FORMAT < IX f ' TYPE OF STRESS IS 'fA4> 

FORMAT ( IXf 12 f 2X f3< E16 ♦ 8f 2X > > 
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055 FORMAT < IX? ' GF 1 ' .5Xr ' SICMAXUJP) ' rHXf ' SIGMAY < GR > ' rQX, 
1 ' S IGMAXY < GF‘ ) ' > 

860 FORMAT ( IX > ' GF" r 5Xf 'EPSLNX(GP) ' rtJXr ' EF'SLNY ( GP > ' r8X» 
# ' EF’SLNX Y < GF’ ) ' ) 

870 FORMAT OX , 'ACCELM KF TURNED UITH ERROR ' r i 4) 

371 FORMAT < IX » ' ACCFES RETURNED WITH ERROR ',14) 

875 FORHAT< IX, ' ACCLCS RETURNED WITH ERROR ' ,M> 

876 FORMAT ( IX, 'tU2DPU RE TURNED WITH ERROR ',14) 

87e FORMAT < IX f ' ACCEI.C RETURNED WITH ERROR ',)4) 

1004 FORMAT (1-1) 

11001 FORMAT (A-l) 

C 

RETURN 

END 



SUBROUTINE ST1 601 ( NT f I r L * F’SIBTB > 

CP**** * ******************* ********************** *********** 

CP* 

CP* ST1601: DESIGN SENSITIVITY VECTOR FOR A BENDING FL AT F 
CP* 

CP************************************************************ 

CP* 

cp* description: 
cp* 


CP* . ' ST 1601 / COMPUTES I ME DESIGN SENSITIVITY VECTOR FOR 
CP* A TRIANGULAR BENDING PLATE ELEMENT# 

CP* 


CP************************************************************ 


CP* 

CP* 

CP* 

CP* 

CP* 


NT COUNTER FOR FINITE DIFFERENCE 

I EXTERNAL ELEMENT NO. BEING PROCESSED 

L EXTERNAL LOAD CASE NOS, 

PSIBTB STRESS CONSTRAINT FOR BENDING PLATE 


CP* 


CP************************************************************* 
INCLUDE ' C AEGSDR . INC 3 IMPLIC.SPC' 

INCLUDE ' C AEGSDR » INC J AOCIPN. MON ' 

INCLUDE ' EAEGSDR . INC 3 CNTL . NON ' 

INCLUDE 7 C AEGSDR . INC3 ELFDES ♦ MON ' 

INCLUDE 'C AEGSDR. 1NCJ SVECTR ♦ MON ' 

COHMON/LCSDES/DLCS( 90 ) 

C 

EQUIVALENCE (NDAT(9S) »IDBL) 


C 

DIMENSION X<3> iY<3)rZ<3> fSXG<;*»3> t l£PN(3r3> fL.<2) f ILCN<2> 9 
* SE ( 500) r BUF (100) FpSIRf B(5u0> fSBUF(50> f RSI BF*£ < 5u0 > 

C 

DATA KT/3/»lREF/l/fMPI/l/f ITYPE/1/ 


C 

C 

SE ( I ) « O.ODO 
PSIBTB<I) = O.ODO 
PSIBPE(I) « O.ODO 


C 


20 DO 50 J J- 1 r 2 


C 

C GET INTERNAL LOAD CASE NUMBER 
C 

CALL ACCLCS< .1 r 3 PNLCSr I DDL f 1 f 0 f I ERR > 

I F < I ERR ♦ NE « 0 ) GO TO 805 

CALL ACCLCS (2f IF’NLLSrL< JJ) 1 2$ DLCS r IERR ) 

IF(IERR.NE.O) GO TO 805 

ILCN(JJ) * DLCS( 21 ) 

C 

C GET PROPERTIES 

C 

CALL ACCEPR ( 2 f lF'NEPR fIPTABfOf BUF f LFN f X F'RR ) 

IF ( I ERR ♦ NE * 0 ) GO TO 809 

PB(NT>=BUF<25) 

IF(NT.GT.l) GO TO 100 
50 CONTINUE 


C 

C GET X AND Y FOR JACOBIAN EVALUATION 


C 

100 CALL ACCELC ( 2f IFNlfLC v KIN f f I REF f BUF fLFNBfIERR) 

IF(IERR.NE.O) GO TO 007 
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M ® 1 

HO 110 J=1 f LENB r 3 
K = J + l 
LL = Ji2 
X ( M ) = BUF < J) 

Y ( M ) = BUMK) 

Z<H> = BUFU.L) 

M * h + 1 

110 CONTINUE 

i; 

C GET STRESSES FROM ORIGINAL LOAD CASE 
C 

CALL SSMD1 6 < X * Y * I LCN (1) * SIG ) 

IF (NT * GT « 1 ) CO TO 120 

C 

C GET STRAINS FROM ADJOINT LOAD 

C 

CALL SNMDl 6 < X * Y * I LCN < 2 ) * EF‘N > 

C 

120 AREA = EUTRIA(XrY) 

XMP = 1, DO/ ARE A 

t; 

C START INTEGRATION LOOP 

C 

no 400 IT=1 *3 

c 

VMS * DSQRT«SIG(ITf 1 > **2+SIG 0 T r 2 > *¥2-SIG < I T * 1 >*SIG< IT* 2) 

* +3*SIG<ITf 3)**2> 

TMAX = DSGRT( < «5*(SlG(lTf l)-SIG(ITr2) ) ) **2+SIG CH , 3 ) **2 ) 
IF(NT.GT.l) CO TO 345 
IF <I.Nt.lt:E(NC)) GO TO 320 
IF (IST.EQ.l) GO TO 300 
IF < 1ST . GT ♦ 2) GO TO 320 
C 

C CALCULATE VON MISES STRESS SENSITIVITY TERM 
C 

SE ( I > =SE < I ) + ( AREA/ 3 ♦ HO > *VMS*XHP/P» (Hi) 

GO TO 320 

C 

C CALCULATE PRINCIPAL STRESS SENSITIVITY TERM 

C 

300 SE ( I >=SE(I H < AREA/3. DO) *< .500* (SIC CM ? 1 ) iSIG ( I T * 2 > +2 

* ♦TMAX ) )*XHP/FB(ND 
C 

C CALCULATE THE SENSITIVITY VECTOR 

C 

320 no 340 0-1*3 

SE(I> - SE(I) - ( AREA/3 • DO >*SIG( IT* J>*EPN<n rJ> 

340 CONTINUE 

BPSITB(I) = SE ( I ) 

C 

345 IF < I * NE ♦ I CF( NC ) ) GO TO 400 

C 

C* CALCULATE PSI(B) - INTEGRAL OF STRESS FUNCTION G FOR ELEMENT 


IF ( 1ST • EQ • 1 ) GO TO 360 

PSIBTB(I) = PSIBTBCI) + (AREA/3 • DO) *VMS*XHP 
GO TO 400 

360 PSIBTB(I) = PSIBTBd >+< AREA/3 . »0> * < .5*<SIG< ITr 1 >+SIG( IT*2> > 

♦ +TMAX>*XMP 

400 CONTINUE 
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c 

c 

GO TO 820 

C WRITE ERROR MESSAGF5 Hi 1 HE SCREEN 
C 

000 PRINT 870 * I ERR 

GO TO 820 

801 PRINT 871* I ERR 

GO TO 820 

805 PRINT 8/5* I ERR 

GO TO 820 

807 PRINT 878* I ERR 

GO TO 820 

80? PRINT 876* I ERR 

GO TO 820 

C 


C 

820 

C 

C 

850 

851 

854 

855 

860 

862 

870 

871 

875 

876 
878 
1004 
2001 
C 


CONTINUE 


FORMAT C IX* ' JOIN I LOAD IS APPLIED AT ELEMENT ' * ) 4 ) 
FORMAT < IX* ' *#*TYPE OF STRESS IS '*A4) 


FORMAT (IX* ]2*2X*J(L16«S*2X)) 


FORMAT (IX* 'GP' * 5X* ' S I UMAX < GP > ' * 8X *' SIGMA Y < GP ) 
*'SIGMAXY<GP) ' ) 

FORMAT< IX* 'CP' * GX* ' FPSLNX ( GP ) ' *8X* 'EPSLNY(CP) 
*'EPSLNXY<GP> ') 


FORMAT (IX*' ELEMEN I '*14) 
FORMAT ( IX* 'ACCELM RETURN!’ u 
FORMAT (IX? ' ACCFES RETURNED 
FORMATdX* 'ACCLCS RETURNED 
FORMATdX* ' ACCEF’R RETURNED 
FORMAT (IX* ' ACCELC RETURNED 
FORMAT ( 1 4 > 

FORMAT <A4) 


WITH 

ERROR 

'* *14) 

W1 \ H 

ERROR 

'*14) 

WITH 

Error 

' *14) 

WITH 

ERROR 

' *14) 

WITH 

ERROR 

' *14) 


* SX * 

* 8X * 


RETURN 


END 
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c 

c 

c 


ELEMENT ATTRIBUTES COMMON 


COHMON/ELEDES/IEfKSO) f 1NTNNC32) 
EQUIVALENCE (IED(l)rl TVP ) 

* ( IED( 4 ) 9 NUOF ) 

* (IED (7) fILUMP) 

* < IED(9) fNRPT) 

* <IEH(11) rKIND 

* <IED(14) f IESM) 

* (IED(17)fIS0P> 

* ( IED (20) f IMA IF') 

* ( IED (23) f NUfiKEI ) 

* ( IED < 28 > r ISPTAB) 

* (IED(31)f NSVAI . ) 


9 1 EXTNN ( 32 ) f 1 1>LD ( A ) , I AF ( 3? ) 


(IED (2) f 
<IED< 5)t 
( IED (8 > f 
< IED(IO) 
(IED( 12) 
( I ED ( 15 ) 
<IED(18) 
(IED (21) 
( I ED (26) 
( I ED (29) 
( I E D ( 3 2 ) 


) STYP > 

MAXDUF) 

iAcru> 

F Kb X r ) 
fKNLX 1 ) 
/N.IND 
fIDOUN) 
f (MATO 
f NhA 1 > 
fDL FA) 

vie r r? \ 


F ( I E D ( 3 ) t NlJNh'E ) F 
f UED(6> 9 NUh.LL ) f 


( IED( 13) fKI. INI ) f 
<IED(16> F JrlCOPI ) F 
( I E D < 1?) f NHL ) f 
( IED (2 2) t hh’EL ) 9 
< I E D ( 2 7 ) fIPTA*) • 

<IFD<43> rh&Iti ) 


COMMON/S VECTR/DPS IT < SOO ) f DPS I D ( SOO ) f DPS IH ( 500 ) f DPSITB ( 500 > f 

* Th( 2 ) f DW(2> f DH(2) fPB(2) f ICT f 1SAC f LCS f Nl C f NC f 

* 1ST f ICE ( 200) 


C 

C *** CURRENT IPNT POINTER STORAGE FOR ACCESS ROUTINE CAMS 
t: 

COMMON/ACC IPN/ IPHKLM ( 2 ) f I PNNOD ( 2 > fIPNNlC(2) f IPN thh<7) *IPNLMM(2 ) f 

* IPNEDK ( 2 ) f I PNELC ( 2 ) f 3 F'NhNb ( 2 > f 3 PNBER f i PNM AT f 

* I PNNCF ( 2 ) f I PNAFL ( 3 ) f 1 PNEEN ( 2 ) f 1 PNAND ( 2 > f \ PN t KR ( 2 ) f 

* IPNNSFC2) fIPNBMF ( 2 ) f 3 PNFES ( 2 ) fIPNLCS(2) f IPNSTK f 

X IPNRES fIPNASD f IPNCLPC 2 > f IPNEPR fJPNMPk f 

It IPNLMI (2) f I PNBLK f IPNCNLi(2) f IPNNDF ( 2 ) f 1PNSGM ( 2 ) f 

* IPNNHF fIPNTEM » XPNISS (2) f IPNSTH 


c 

C*** GLODAL CONTROL PARAMETERS 
C 

COMMON /CNI L/ NETY<200) fNDAT<300) f NL INE f NWID4 f Null D1 f NAUX f 
X IHEADR < 1 32) f IAUX(33f5> 



